NAS A-RP- 1153 19860013108 


NASA 

Reference 

Publication 

1153 


April 1986 


A Methodology for Airplane 
Parameter Estimation and 
Confidence Interval 
Determination in Nonlinear 
Estimation Problems 


Patrick C. Murphy 




j 


LAf.'CLLY n-"S£. 
UPPARY 
HAMPTON, 



CfJ 





NASA 

Reference 

Publication 

1153 


1986 


A Methodology for Airplane 
Parameter Estimation and 
Confidence Interval 
Determination in Nonlinear 
Estimation Problems 


Patrick C. Murphy 

Langley Research Center 
Hampton, Virginia 


NASA 

National Aeronautics 
and Space Administration 


Scientific and Technical 
Information Branch 



SUMMARY 


A methodology is presented for airplane parameter estimation and confidence interval 
determination in nonlinear estimation problems. In addition, as part of the methodology, 
an efficient scheme to determine aerodynamic model structure is suggested and briefly 
described. The algorithms described provide a unified approach to solving nonlinear airplane 
identification problems. The nonlinear estimation problems of interest to this study are 
further complicated by nonlinear system dynamics and in particular nonlinear aerodynamic 
models. 

An algorithm for maximum likelihood (ML) estimation is developed with an efficient 
method for approximating the sensitivities. The algorithm is applicable to most parameter 
estimation problems and is particularly suited for nonlinear, multivariable, dynamic systems. 
The ML algorithm relies on a new optimization method closely related to a modified Newton- 
Raphson (MNR) technique; the new method is referred to as modified Newton- Raphson with 
estimated sensitivities (MNRES). 

MNRES determines sensitivities by using slope information from local surface approxi- 
mations of each output variable in the parameter space. The fitted surface allows sensitivity 
information to be updated at each iteration with a significant reduction in computational 
effort. With MNRES, the sensitivities can be determined with less computational effort than 
with either a finite-difference method or integration of the analytically determined sensitivity 
equations. The type of surface (for example, nth-order polynomial or spline) and the method 
of fitting the surface (for example, least squares or solution of simultaneous equations) are 
chosen by the user to suit the particular need. MNRES eliminates the need to derive sen- 
sitivity equations for each new model, thus eliminating algorithm reformulation with each 
new model and providing flexibility to use model equations in any convenient format. 

Two surface-fitting methods are discussed and demonstrated, while other possibilities 
are indicated. MNRES is compared with other commonly used optimization methods, a 
search method called the flexible polyhedron search (FPS) and a gradient method called 
the modified Newton-Raphson method. Several sample problems are solved to compare the 
techniques. Simple linear systems are used at first, and then nonlinear aircraft estimation 
problems are solved by using both real and simulated data. MNRES is found to be equally 
accurate and substantially faster than the commonly used techniques. The reduction in 
computational effort provided by MNRES depends on several factors: the choice of surface- 
fitting method, the number of unknown parameters, data quality, accuracy of the sensitivity 
calculations, and, particularly, the degree of nonlinearity of the cost function. 

A search technique for determining the confidence limits of ML parameter estimates is 
applied to nonlinear estimation problems for airplanes. The confidence intervals obtained 
by the search are compared with Cramer-Rao (CR) bounds at the same confidence level. It 
is observed that the degree of nonlinearity of the cost function is an important factor in the 
relationship between CR bounds and the error bounds determined by the search technique. 
The CR bounds were found to be close to the bounds determined by the search when the 
degree of nonlinearity was small. The CR bounds were 3 to 8 times too conservative (too 
small) when the nonlinearity was significant. Beale’s measure of nonlinearity is developed 
in this study for airplane identification problems; it is used to empirically correct confidence 
levels for the parameter confidence limits. The primary utility of the measure, however, was 
found to be in predicting the degree of agreement between CR bounds and search estimates. 
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I. INTRODUCTION 

Problems in dynamics may be divided into three 
categories. By considering a general dynamical sys- 
tem f with input U and output Y the categories can 
be defined as (1) the classical problem where U and f 
are given and Y is to be determined, (2) the controls 
problem where f and the desired Y are given and U 
is to be determined, and (3) the identification prob- 
lem where U and Y have been measured and f is to 
be modeled. 

The theory for system identification provides a 
way for modeling an unknown system on the basis 
of input and output information. Identification the- 
ory incorporates a priori knowledge of the dynamic 
processes and stochastic processes involved; thus, the 
identification problem is not usually characterized as 
a “black box” problem. In fact, system identifica- 
tion problems are usually characterized, as in refer- 
ence 1, by three factors: (1) class of models, (2) class 
of inputs, and (3) a criterion for state and parame- 
ter estimation. The models and inputs may be de- 
terministic or stochastic and the criterion (cost func- 
tion) may be based on statistical theory or numerical 
considerations. 

Implementation of identification theory usually 
follows four basic stages. The first stage requires 
the design of an experiment, that is, specification 
of the identification objectives, statement of system 
configuration and conditions, and selection of an in- 
put form. Determining an optimal input for iden- 
tification can be critical; all the modes of a system 
must be excited in order to identify the system cor- 
rectly and completely. The second stage is model 
structure determination (a more comprehensive term 
is model characterization). The model is assumed 
to be linear or nonlinear, time varying or time in- 
variant, with or without process noise, with or with- 
out measurement noise, etc. The unknown parame- 
ters in the model may include system parameters as 
well as initial conditions, bias terms, and character- 
istics of measurement and process noise. The third 
stage involves parameter and state estimation, which 
provides mean values and standard error estimates. 
These are obtained by finding an extremum of some 
optimality criterion. State estimation can be better 
characterized as a filtering problem; a Kalman filter 
is commonly used. The fourth stage is verification, 
accomplished by comparing estimates from experi- 
mental data sets and other estimation techniques. In 
addition, other sources provide comparisons; for air- 
planes, both wind tunnel and theoretical predictions 
are used. Verification is also accomplished through 
sensitivity analysis and through analysis of residuals 
and model predictive capabilities. 


System identification theory has become impor- 
tant to aircraft technology for several reasons. It pro- 
vides an alternate approach to determining aircraft 
characteristics (parameters). Comparing results with 
other techniques is always good scientific practice. 
Purely theoretical approaches or purely experimen- 
tal approaches (wind tunnels) have in many instances 
failed to accurately predict prototype characteristics. 
By offering an opportunity to observe actual vehicle 
performance, flight testing results in better calibra- 
tion and understanding of wind tunnel results and 
more accurate modeling for ground-based simulators. 

The development of aircraft parameter estimation 
paralleled developments in estimation and system 
theory. Early flight test studies centered on steady- 
state maneuvers and free oscillations. These studies 
were time consuming and provided limited informa- 
tion. The main interest was to obtain basic aero- 
dynamic parameters, termed stability and control 
derivatives, from linear dynamic models combined 
with linear aerodynamic models. In the early 1950’s 
Greenberg and Shinbrot developed a least squares 
approach to analyzing simple transient maneuvers 
(refs. 2, 3, and 4). However, without computers the 
simplest flight test problem with only four unknown 
parameters took 24 hours to analyze (ref. 5). A ma- 
jor development for aircraft parameter estimation oc- 
curred in the mid 1960’s. Large-capacity, high-speed 
digital computers and highly automated data acqui- 
sition systems were introduced. In 1968, when Lar- 
son (ref. 6) applied the method of quasi-linearization 
and Taylor and Iliff (ref. 7) introduced the modified 
Newton-Raphson method, a new stimulus was given 
to parameter estimation. Other contributions came 
in the early 1970’s from Mehra (ref. 8), Stepner and 
Mehra (ref. 9), and Rault (ref. 10). 

During the past decade application of estimation 
theory to nonlinear systems has become an increasing 
concern, stimulated primarily by aerospace applica- 
tions. Today, nonlinear dynamics is commonly incor- 
porated with linear aerodynamic models in flight test 
data analysis. The techniques are well established for 
flight regimes where the aircraft aerodynamic model 
can be expressed as a linear function of states and 
control inputs. However, modeling the combination 
of nonlinear dynamics and nonlinear aerodynamics 
and estimating the parameters associated with that 
model present many difficulties. The need to identify 
the best mathematical representation (model struc- 
ture) and estimate the associated parameters for non- 
linear flight regimes has motivated further develop- 
ment of identification and estimation techniques. 

A new approach to airplane parameter estimation 
and confidence interval determination is offered in 
this study as a contribution toward building a more 



general and unified airplane identification methodol- 
ogy. The more general methodology starts with the 
work done in reference 11, which suggests a useful 
technique for model structure determination when 
nonlinear aerodynamic effects are present. The sug- 
gested technique uses a modified stepwise regression 
(MSR) along with several testing criteria to deter- 
mine a parsimonious, yet adequate, model. The lim- 
itation of this technique (as with any least squares 
method) is that the estimates are asymptotically bi- 
ased and variance estimates are based on simplifying 
assumptions that are valid only for the “classical” lin- 
ear regression. This limitation can be skirted by ap- 
plying the commonly used maximum likelihood (ML) 
technique, using the model structure determined by 
the regression and the regression estimates as an ini- 
tial guess. The ML approach has much more favor- 
able asymptotic properties (ref. 12), and it provides 
estimates of the Cramer-Rao (CR) bounds for the 
parameter variance. 

There is a computational cost, however, for the 
more favorable asymptotic properties of the ML tech- 
nique. Dynamic systems, such as aircraft, require 
substantial computational effort at each step of the 
optimization process. At each step the equations of 
motion must be integrated to obtain time histories of 
each state and output variable. In addition, most ML 
algorithms use a modified Newton-Raphson (MNR) 
optimization scheme, which requires integration of 
sensitivity equations. This accounts for most of the 
computational effort, since the number of state and 
sensitivity equations to be integrated at each itera- 
tion is equal to the number of states plus the product 
of the number of states and the number of unknown 
parameters. Several states and 20 to 30 unknown pa- 
rameters are not uncommon for one flight condition. 
If a model is desired throughout the entire flight enve- 
lope, the computational requirements become over- 
whelming, particularly when analysis of various flight 
conditions requires more than one candidate model. 
A very efficient ML estimation algorithm is desirable 
to reduce the computational requirements for pro- 
cessing a large number of parameters and candidate 
models. 

Besides the high computational cost associated 
with the ML/MNR algorithm, an additional diffi- 
culty is that it requires the user to have prior knowl- 
edge of the model structure to formulate the sensitiv- 
ity equations and, thus, to formulate the algorithm. 
This can be very burdensome when modeling aircraft 
in nonlinear flight regimes, since model structure may 
change significantly from one flight condition to an- 
other. Therefore, an algorithm that is independent 
of sensitivity equations is very advantageous. 


Reducing computational requirements of the ML 
method requires careful examination of the opti- 
mization methods used in the algorithm. Although 
nonlinear, unconstrained optimization problems have 
been studied quite extensively (ref. 13), little has 
been done to improve the optimization techniques as 
they apply to aircraft estimation problems. Gupta 
and Mehra (ref. 14) considered the numerical aspects 
of computing ML estimates for linear dynamic sys- 
tems in state-vector form and methods for speeding 
up convergence. Bowles and Straeter (ref. 15) con- 
sidered computational aspects for several methods 
used in aircraft identification. Trankle, Vincent, and 
Franklin (ref. 16) considered the difficulties associ- 
ated with use of a nonlinear dynamic model in ML 
parameter estimation and parameter covariance es- 
timation; sensitivity calculation methods were also 
considered. More recently, Trankle, Vincent, and 
Franklin (ref. 17) considered the overall methodol- 
ogy of system identification for nonlinear aerody- 
namic models including computational aspects of the 
problem. In reference 18, a nonlinear least squares 
algorithm is developed which uses a linear-surface 
approximation of a scalar response variable to elim- 
inate derivative calculations altogether. The algo- 
rithm is tested on problems that do not involve 
dynamic systems. Presented in the current study 
is a significantly improved maximum likelihood al- 
gorithm, initially developed in reference 19, which 
relies on an optimization scheme referred to as a mod- 
ified Newton-Raphson method with estimated sensi- 
tivities (MNRES). A surface approximation is also 
used in MNRES; however, it is treated differently by 
developing an algorithm that retains derivative infor- 
mation in a Newton-Raphson method for multivari- 
able, dynamic systems. This is done to provide direc- 
tional information for the convergence process and to 
provide covariance information. With the MNRES 
approach, sensitivity equations are eliminated and 
computational demand is significantly reduced. 

Another difficulty in using the ML technique is 
that the CR inequality provides only a lower bound 
measure of precision for an unbiased estimator. It 
is known from practical application of ML that this 
lower bound can differ from the variance obtained, 
for example, by repeated measurements (ref. 20). At- 
tempts have been made, therefore, either to mod- 
ify the CR bounds by considering a band-limited 
measurement noise (refs. 21 and 22) or to estimate 
the parameter variance directly from measured data 
(ref. 10). Advances in statistical methods also came 
about with the availability of high-speed computers. 
In 1960, Beale (ref. 23) considered the problem in 
nonlinear estimation of determining the approximate 
parameter confidence regions using likelihood ratios. 
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In addition, a measure of nonlinearity was developed 
to assess the quality of the approximation. Surpris- 
ingly, Beale’s work has had very little application 
since it was published (ref. 24). In 1979, Mereau and 
Prevost (ref. 25) used the likelihood ratio approach 
to determine confidence regions for aircraft systems. 
In 1980, Mereau and Raymond (ref. 26) developed a 
search procedure to find the “iso-distances” defining 
the confidence regions. 

The goal of the current study is to provide a 
unified methodology for solving nonlinear airplane 
identification problems. Improved techniques for es- 
timating parameters and their confidence limits in 
nonlinear, multivariable, dynamic systems, in partic- 
ular, aircraft systems, are presented. The improved 
techniques provide (1) increased efficiency for the 
estimation process, (2) elimination of the need for 
a priori knowledge of sensitivity equations, (3) more 
accurate assessment of the parameter error bounds 
than those obtained using CR bounds, and (4) an 
adaptation of Beale’s approach to the airplane esti- 
mation problem. 


The development of this work begins in section II 
with a description of the airplane model and the re- 
gression method used to determine it. Section III 
describes the parameter estimation techniques used 
and their statistical properties. The primary estima- 
tion method used in this study to determine airplane 
parameters is maximum likelihood. Linear regression 
methods are also presented, since special forms are 
developed for use in MNRES and because they are 
used in the model structure determination scheme 
of section II. In section IV the MNRES method is 
presented with a discussion of its various forms and 
properties. Also presented briefly are some com- 
monly used optimization methods, which are com- 
pared with MNRES. Section V develops the theory 
for confidence interval estimation and the adaptation 
of Beale’s work to the airplane problem. Finally, sec- 
tion VI presents the results and discussion for the 
application of these methods to simulated and real 
data. 
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II. MODEL CHARACTERIZATION 

Model characterization, as discussed in the intro- 
duction, establishes the known or assumed charac- 
teristics of the model to be used in the identification 
process. Since system identification usually does not 
involve a black box problem, where nothing is known 
about the model in advance, qualitative statements 
describing the class of model, optimal inputs, and 
statistical properties of the measurements are nor- 
mally provided. 

Generally, the more information available to char- 
acterize the model, the greater the likelihood of suc- 
cessful identification. Of course, attempting a com- 
plete representation of a dynamic system, such as 
an airplane, is extremely difficult, if not impossible. 
Actually, a complete model is unnecessary. The ob- 
jective in identification is to select the simplest model 
that allows proper determination of the desired un- 
known parameters from measured data. The princi- 
ple of parsimony is usually applied. This principle 
states that given a choice of two models having equal 
residual variances, choose the model with fewer pa- 
rameters. Therefore, the objective in identification is 
to choose a parsimonious, yet adequate, model. Very 
complex models may be justified to obtain an accu- 
rate description of the system motion, but it is clearly 
detrimental to the estimation process. If the infor- 
mation content in the measured data is very limited 
or if too many parameters are required, the estima- 
tion algorithm may provide inaccurate estimates or 
it may fail. 

The models considered in this study represent 
dynamic systems, characterized by having derivatives 
with respect to time included in the model in addition 
to the dependent and independent variables. One of 
the possible general forms for these systems is 

X s = /(X s ,U,0, t) X s (0) = X s0 (2.1) 

Y = /i(X s ,U,0,f) (2.2) 

where X s is a vector of state variables, Y is a vector 
of output variables, U is a vector of input variables, 
and 0 is a vector of unknown parameters. The time 
variable t may or may not appear explicitly. This 
form is not as restrictive as it first appears; many 
problems can be cast into the matrix differential form 
above. 

Several difficulties arise when estimation tech- 
niques are applied to dynamic systems. A major dif- 
ficulty is the significant computational demand asso- 
ciated with solving matrix differential equations. In 
an estimation algorithm these equations are solved 
repetitively. A difficulty can also arise when inte- 
grating the equations of motion because the bound- 



Figure 1. Airplane body-axis system with forces and 
moments. 


ary conditions or initial conditions are not always 
known exactly. Therefore, in many estimation prob- 
lems the initial conditions are treated as unknown 
parameters. Another difficulty is that the solution 
to the differential equations can vary depending on 
sometimes very small changes in the unknown pa- 
rameters. For example, a first-degree system is sta- 
ble or unstable depending only on the sign of the 
damping term. Nonlinear systems can amplify this 
type of problem. The success of the estimation may 
depend on the initial guesses for the parameters be- 
cause failure may occur when a parameter is outside 
a stability boundary. Unfortunately, obtaining sta- 
bility boundaries is really practical only for linear, 
time-invariant systems. Finally, numerical difficul- 
ties with truncation and rounding errors are always 
present when numerical differentiation and integra- 
tion are performed. 

A. Airplane Equations of Motion 

The particular dynamic system of interest to this 
study is the airplane, modeled by equations in the 
general form 

X s = /(X s , U, 0) X s (0) = X s0 (2.3) 
Y = h{X s ,V,0) (2.4) 

The equations of motion used are referred to a body- 
axis system (fig. 1). The equations were developed 
with the following assumptions: 

1. The airplane is a rigid body. 

2. The effect of spinning rotors is negligible. 

3. The airplane has a plane of symmetry in the 
XZ- plane. 

4. There are no external disturbances to the air- 
plane. 
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5. Thrust is accounted for as part of Cz and Cx 
where 

Cx = Ct cos a j 1 + C E sin a — C E cos a 
Cz = —Cq 1 sin a j 1 — Cl, cos a — C E sin a 

The resulting nine equations represent a six- 
degree-of-freedom, coupled, nonlinear system where 
the kinematic relations are expressed in terms of di- 
rection cosines. These equations of motion are given 


as follows: 




u = 

-qw 4- rv + qtxz 

:+^C X 

m 

(2.5) 

V = 

—tu + pw + gt yi 

, + ^FL c y 

m 

(2.6) 

U) — 

-pv + qu + g£ zz 

+ ^C Z 
m 

(2.7) 

P = 

F i + 

r r __ r2 1 ^ 

lx*z — izx 

Izx r, 

Ixlz-I'L 

(2.8) 

Q = 

Izx t 2 2\ , 

-") + 

pU -I x + qS w -C Cm 
1 y 1 y 

(2.9) 

r = 

I I 1 " 12 F * + 

lX*Z *zx 

Izx F, 
Ixlz-Ilx 

(2.10) 

ixz = 

r£yz - ql zz 


(2.11) 

II 

— fizz + p£zz 


(2.12) 

&zz — 

qtxz — p£yz 


(2.13) 

where 




Fi 

= ( h ~ Iz)qr + 

IzxPQ + QSyjbwCi 

(2.14) 

f 2 

= (lx ~ Iy)pq ~ 

IzxQT + qSwbwGn 

(2.15) 


The nondimensional aerodynamic forces and mo- 
ments, Cx, Cy, Cz, Ci, C m , and C n (shown in 
fig. 1), are usually approximated by a Taylor series 
expansion around steady trimmed flight conditions 
or by polynomial splines (see ref. 27). The form of 
the aerodynamic model equations is 

y{t) = 00 + 0l x l + @2 X 2 H + 0n p -l x n p -l (2-16) 

or in vector form 

= (2.17) 

where y(t) or j/j represents one of the nondimensional 
aerodynamic forces or moments at time t or at the 
ith data point. The stability and control derivatives 
are represented by 6\ to 0 np - i, and the trim forces or 


moments corresponding to an initial trimmed flight 
condition are represented by 0q. The x\ to x np _i 
represent any functions of the state and control vari- 
ables chosen for the model. The row vector X t - is 
given as 

Xj- = [l xi x 2 x np _i] (2.18) 

In general, the form of the aerodynamic equation is 
unknown; however, for estimation it must be pos- 
tulated. The form may vary significantly from one 
flight condition to another. 

The following output equations, used in this 
study, reflect the measurements commonly available 
from flight tests: 


V = Vu 2 + v 2 + w 2 

(2.19) 


(2.20) 

° = tai1 " 1 o 

(2.21) 

0 = sin _1 (-4 2 ) 

(2.22) 


(2.23) 

a x = + qw — rv — gt xz ) 

(2.24) 

a V = “(*> + ru-pw- gtyz) 

(2.25) 

1 

a z = ~(w + pv — qu — gt zz ) 

(2.26) 

The airplane identification problem can be made 
more tractable by treating longitudinal and lateral 
cases separately. This is accomplished by providing 
the required lateral information to the longitudinal 
case (or the required longitudinal information to the 
lateral case) in the form of measured input variables. 
This has been used successfully in many other stud- 
ies, for example, reference 20. Thus, the states, out- 
puts, and inputs for the two cases are given as follows: 

For the longitudinal case, 


X s = [u w q £ xz £yz £zz\T 

(2.27) 

Y = [V a q 0 a x a z ] T 

(2.28) 

U = [8 e j3 E v E p E r E <f> E ] T 

(2.29) 
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For the lateral case, 

X = [v p r i xz ^yz £ zz\ (2.30) 

Y = [/3 pr cj> a y ] T (2.31) 

U = [<5 0 8 r ue we Qe Qe a E] T (2.32) 

where the subscript E indicates a measured quantity. 

B. Model Structure Determination 

The goal of model structure determination is to 
determine an analytical representation of the sys- 
tem that can be classified as adequate. An adequate 
model is one which sufficiently fits the data, allows 
successful estimation of the parameters, and has good 
prediction capabilities. In aeronautical applications 
the form of the rigid body equations of motion is 
known. The primary uncertainty, with regard to 
model structure, is in the aerodynamic model equa- 
tions (eq. (2.16)). One of the successful methods for 
determining the model structure of these equations 
from measured data is based on stepwise regression. 

In the stepwise regression approach, after postu- 
lating the aerodynamic model equation, significant 
terms among the candidate variables are determined 
and corresponding parameters are estimated. The 
variable chosen for entry into the regression equation 
is the one that has the largest correlation with y after 
adjusting for the effect on y of the variables already 


selected. The parameters are estimated by the least 
squares technique. At every step of the regression, 
a new variable entering the model and the variables 
incorporated into the model in previous stages are 
reexamined. Any variable providing a nonsignificant 
contribution (due to correlation with more recently 
added terms) is removed from the model. The pro- 
cess of selecting and checking variables continues un- 
til no more variables are admitted to the equation 
and no more are rejected. Experience shows, how- 
ever, that the model based only on the significance of 
individual parameters in model equation (2.16) can 
still include too many terms and, therefore, may have 
poor prediction capabilities. Several criteria for the 
selection of an adequate model are introduced in ref- 
erence 11 and the details of the whole procedure are 
explained in references 11 and 28. Two procedures 
utilizing the algorithm in reference 11 offer certain 
advantages by capitalizing on the use of splines for 
modeling (ref. 27) and large amplitude maneuvers for 
efficient data analysis (ref. 29). 

The stepwise regression procedure in reference 1 1 
provides a very efficient method of determining aero- 
dynamic model structure and initial parameter esti- 
mates. This makes the computationally demanding 
problem of modeling the entire flight envelope more 
tractable. Once the aerodynamic model structure is 
determined, a more advanced method is needed to 
improve the biased parameter estimates obtained in 
the regression. More advanced methods are discussed 
in the next chapter. 
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III. PARAMETER ESTIMATION 

In this study maximum likelihood (ML) and lin- 
ear regression (LR) techniques are used to estimate 
parameters. ML is used to estimate both airplane 
parameters and their standard errors (Cramer-Rao 
(CR) lower bounds) from flight data. The ML al- 
gorithm is used with various optimization schemes, 
which are described in section IV. LR is used for 
three different applications in this study: (1) esti- 
mating aerodynamic model structure, (2) estimating 
airplane parameters (starting values for ML), and 
(3) estimating sensitivities in MNRES. The first and 
second applications of LR were accomplished using 
stepwise regression as described in the previous sec- 
tion. The third application was accomplished using 
an algorithm developed in this study. 

A. Linear Regression 

Linear regression analysis is a part of statistical 
theory that generally deals with the determination 
of relationships between response and predictor vari- 
ables. One application of LR theory is curve fitting 
or surface fitting. In this application, the predic- 
tor variables (independent variables) are assumed to 
be deterministic and known without error; response 
variables (dependent variables) may have error. A 
numerical method commonly used in curve fitting 
to compute empirical coefficients is the method of 
least squares (LS). In this method, the same model 
form as equation (2.16) can be used to fit the curve 
or surface. The solution for the unknown parame- 
ters or coefficients is found by minimizing the sum of 
squares of the error between known data points and 
computed data points determined by the model. The 
LS method is valid only for linear problems, that is, 
problems where the unknown parameters occur lin- 
early in the model regardless of whether the model 
structure itself is linear or nonlinear. A LS problem 
can be solved in a batch mode or recursive mode and 
both modes have application for determining the sen- 
sitivities in MNRES. A comprehensive discussion of 
regression analysis is given in reference 28. 

1. Batch Processing 

Batch processing of data in the LS method is 
probably the most commonly used approach for 
curve-fitting problems. The model form given by 
equation (2.17) can be written as 

yt=X,-S + e f {i = 1, 2, N) (3.1) 

where is the fth value of one response variable; 
Xf is the zth vector of predictor variables; S is 
the vector of unknown coefficients; and e; is the 


equation error at the fth data point. This error may 
contain measurement noise, process noise, and/or 
modeling error. However, no assumptions are made 
about the statistical properties of e { (see ref. 28). In 
application to MNRES, y{ represents one element of 
the output vector Y(0); X 2 - represents the ith set of 
values for the vector of unknown parameters 0; and 
S is the np-dimensional vector of coefficients to be 
computed. If a first-degree n p -polynomial expansion 
is chosen for Xj, then the elements of S are the 
desired sensitivities (slopes). 

Applying the LS criterion, which requires min- 
imization of the mean square error, gives the cost 
function as 

N N 

J(S) = ]Te! = ^(y ; -X i S) 2 (3.2) 

i=i l=i 

and minimization requires that 

- 0 = £ xf (Vi - XjS) (3.3) 

1=1 

Solving for S gives 

( N \ -1 N 

s= xx x . p' 4 ) 

Vi=l 7 i=l 


2. Recursive Processing 

Recursive processing significantly reduces mem- 
ory requirements for the MNRES algorithm. How- 
ever, a specialized form of recursive least squares is 
needed for surface fitting in MNRES. Normally in a 
recursive LS problem the purpose is to update pa- 
rameter estimates based on N data points with some 
new information, so that the updated estimates are 
based on N + 1 data points. In the following deriva- 
tion a LS recursive algorithm is designed specifically 
for the MNRES algorithm. MNRES requires the pa- 
rameters to be updated by including new informa- 
tion while removing old information, so that the es- 
timates are always based on a constant number of 
data points. 

As in the batch mode, the surface fitting is per- 
formed to obtain slope or derivative information. 
Consider the LS problem formulated as 

Y = XS + e (3.5) 

where Y is a vector of n data points on a surface to 
be fit by the model given as XS, S is a vector of n p 
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unknown coefficients (slopes), and X is an n by n p 
matrix defined as 

1 Xu X 12 ‘ ' ' X\ np 

1 £21 x 22 x 2n v 

: : ; : (3-6) 

1 X n l X n2 • • • X nrl p _ 

LS estimation gives a solution as 

S = (X r X)- 1 X r Y (n > n p ) (3.7) 

Now define a recursive relation for the ( k + l)th 
iteration 

Pfc+i — (3.8) 

and the updating equation as 

Pfc+l = (Pfc 1 + a£a k ~ a° k ) * (3-9) 

where the row vector a*, is the new set of 0 to be 
included in X and a£ is the outgoing set of 0, which 
produced the highest value of the cost function. The 
recursive relation for S is 

Sfc+i = S k — Pfc+i [a£ y k - a ' k y k 

+ ( a k a k' a f a °k) S k\ (3-10) 

where y k and y k are the outgoing and incoming scalar 
elements of the vector Y, respectively, at the Arth 
iteration. 

The derivation for equations (3.9) and (3.10) is 
as follows. Define Z to be the common elements of 
X between two iterations. Partitioning X for the 
(k— l)th and /cth iterations results in 


a o -| 

X*-1= * 

l^ki 

(3.11) 

-y a A: 

Xfc = rj 

V L k\ 

(3.12) 

From equations (3.11) and (3.12), 

Xfc_i X fc _! — a£ a£ + Z T k Z k 

(3.13) 

XfcXfc = a^a k + ZfcZ k 

(3.14) 

From equation (3.13), 

z T k z k = x T k _ l x k _ 1 - a ° ka ° k 

(3.15) 


Substituting equation (3.15) into (3.14) gives 

XfcXj. = XI’-jXj.-j — a£ a^ + aja*. (3.16) 
Substituting equation (3.16) into (3.8) gives 

Pfc+i = (Xfc-iXfc-l - a °k a °k + aj[a fc ) -1 (3.17) 

which can also be written as equation (3.9). 

Applying the same development to equation (3.7) 
gives 

Sfc+i = Pfc+i (xLjYfc.! — a k y k + a^y^j 

(3.18) 

and substituting equations (3.7) and (3.8) delayed a 
step into equation (3.18) gives 

Sfc+i = Pfc+l ( p fc ls fc - af y% + a (3.19) 

Expanding equation (3.19) gives 

Sfc+l = Pfc+iPfc - P k+1 af y% + P fc+1 a ^y k 

(3.20) 

Note that 

S k - P k+1 P^ t S k = 0 (3.21) 

and then add equation (3.20) to (3.21) to obtain 

Sfc+i = S k — Pfc + iP + P k+1 P^S k 

~ Pfc+l a fc y k + P k+l a k y k (3.22) 

Combining terms in equation (3.22) gives 

Sfc+i = S k — Pfc+i [(Pj.+i — P fc S k 

+ *fv 0 k-*kVk] (3.23) 

and substituting equation (3.9) into (3.23) yields 
equation (3.10). 

3. Statistical Properties of LS Estimates 

Although the LS technique, a numerical proce- 
dure, is not necessarily based on any statistical for- 
mulation, the LS estimator is often characterized in 
statistical terms, since the estimates can be treated 
as random variables. 

In the general LS problem, both process noise and 
measurement noise occur in the data. The model has 
the form 

Y t =X t 0 + w p (3.24) 
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where Yt is an N by 1 vector of the response variable 
for the N data measurement points, subscript t 
indicating that this is the true value of the variable 
without measurement noise, X( is an N by n p matrix 
of the state and input variables; 0 is the n p by 1 
vector of unknown parameters; and w p is the N by 
1 vector of process noise. The measurements provide 

Y = Y t+Vy (3.25) 

X = X t + v x (3.26) 

where v y and v x are the measurement noise in Y 
and X. The process noise and measurement noise are 
typically assumed to be stationary, zero mean, and 
independent random processes. The solution to this 
LS problem (from the last section) is 

0 = (X r X) _1 X T Y (3.27) 

Premultiplying equation (3.24) by [X r X] -1 X r and 

substituting equations (3.25), (3.26), and (3.27) re- 
sults in 

0 = 0 + (X r X)- 1 X T ( wp + vy — v x 0) (3.28) 

Therefore, the expected value of the estimate error 
has the form 

E{0 -0} = -E{{X T X)- l X T v x }0 (3.29) 

and the covariance matrix is 

CO v(9-0) = E{(X T X)- 1 X T ee T X(X T X)- 1 } 

+ E{CX t X.)~ 1 X t v x 
00 t v i x(x t x)- 1 } (3.30) 

where 

e = w p + vy 

From these equations it can be seen that the LS 
estimator is biased even if the measurement noise, v x 
and v y , and the process noise, w p , are zero mean and 
independent. Only with the additional assumption 
that X is known without error (i.e., v x = 0), as 
might be the case in a curve-fitting problem, will 
the estimates be unbiased. Note that the covariance 
matrix is affected by all the measurement and process 
noise. 

B. Maximum Likelihood 

General parameter estimation for an airplane in- 
volves solving the nonlinear estimation problem in 
the presence of both process and measurement noise 


while modeling the airplane with the coupled, nonlin- 
ear equations of motion. One of the advanced tech- 
niques commonly used for this problem because of its 
superior statistical properties is maximum likelihood 
(ML). 

1. Algorithm Development 

Assume that the outcome of an experiment is N 
observations of the n 0 by 1 output vector Zj for 
i = 1, 2, . . . , N, which depends on the unknown pa- 
rameter vector 0. In general, the unknown param- 
eters are the aerodynamic parameters, initial condi- 
tions, and measurement and process noise statistics. 
Let /(Zi, Z 2 , . . . , Z^r/0) be the conditional proba- 
bility density function for the measurements given 
0. The maximum likelihood estimate is the esti- 
mate for which the outcome of the experiment Z t - 
for i = 1,2,..., TV is most likely to occur; that is, the 
probability density is maximized with respect to 0. 
The ML estimate can be expressed as 

0 5 max 9 /(Z 1 , Z 2 , . . . , Z N /0) = /(Zj , Z 2 , . . . , Z N /0) (3.31) 

Using the property of joint conditional probability 
density functions that 

/(Zj, Z 2 , Z 3 /0) = f(Z 3 /Z 2 ,Z 1 ,9)f(Z 2 /Z u 0)f(Z 1 /0) 

(3.32) 

allows the density function to be written as 

/( Zi, Z 2 , . . . , Z N /0) = /( Z n /Z' n _ v 9) 

x /(^AT — l/Z)y_ 2 , 9) 

x ■ ■ ■ f{z 2 /Zi,0)f(Zi/0) (3.33) 

N 

= Y[f(Z i /z' i _ v 0 ) (3.34) 

i=i 

where 

Zi-i = Z i _ 1 ,Z i _ 2 ,...,Zi (3.35) 

If the Zj measurements are treated as fixed, the den- 
sity function becomes a function of 0 only. This func- 
tion is usually referred to as the likelihood function, 
that is, 

N 

L l (O) = Hf(Z i /Z' i _ l ,0) (3.36) 

i— 1 

Consequently, the problem of finding a maximum 
likelihod estimate becomes the problem of finding the 
0 that maximizes the likelihood function. 

To define the likelihood function, the distribution 
of the measurements given 0 must be defined. If the 
distribution of the measurement and process noise is 
Gaussian, then the distribution of the measurements 
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given 0 is also Gaussian and can be uniquely deter- 
mined by computing the mean and covariance. The 
mean is 


E{Z i /Z' i _ v O} = Y i/i _ l (0) (3.37) 


where Y i/,-_i(0) is the best estimate of the measure- 
ment at time t,- given measurements up to and includ- 
ing the previous point. By definition, the covariance 
is 


cov{Zi/Z •_!,#} = E{(Zi - Y^XZi - Yj / j_ 1 ) 7 '} 

= E{v iV J} = Bj/j.jtJ) (3.38) 


where v ,■ is the vector of residuals. Now the problem 
is to compute the conditional mean Y and 
covariance B t yj_ 1 (0). For systems including process 
noise these values can be obtained with a Kalman 
filter. It has been shown (ref. 8) that for high 
sampling rates (as are commonly used to collect 
flight test data), the residuals tend toward a 
Gaussian distribution. Therefore, the distributions 
for both ly and (Z t -/Z t -_i, 0) are reasonably assumed 
to be Gaussian. In systems without process noise, 
some simplifications are possible. In particular, the 
residual error may be written as 


v i = Z i -Y i 0) (3.39) 


where Y,- is the predicted value of the output vec- 
tor at time t{. Without process noise the conditional 
subscript is not needed since the Kalman filter up- 
date is zero. Also, assuming stationary statistics, the 
mean and covariance of the residuals may be written 
as 


E{ui) = 0 

(3.40) 

E{vivJ } = R 8ij 

(3.41) 


where 8{j is the Kronecker delta. With these as- 
sumptions, the conditional probability density func- 
tion can be written as 


/(Z i /Zj_ 1 ,tf)(= (27r)- n °/ 2 |R|- 1 / 2 exp V t ) 

(3.42) 

Hence, using equation (3.42), the likelihood function 
is 


Li(e) = Hf(z i /z > i=1 ,0) 

i = 1 

= (2n)~ Nn °/ 2 JI |R| -1 / 2 


X exp 


(3.43) 


The negative logarithm of the density function is 
more convenient to use and since the density function 
is non-negative, there is no change to the problem 
except that maximizing the density function equates 
to minimizing the negative logarithm of the density 
function. Thus, the negative log likelihood function 
is more commonly used: 


i W = |D Z i-Y,O r R" 1 (Zi-Y f ) 

L i = 1 

N 

+ — In |R| + Constant (3.44) 

4U 

The unknown R can be estimated by minimization 
of the likelihood function with respect to R. This 
produces 


6 -=^E( Z i-Yi)(Z i -Y i ) T 

i = 1 

= ^J2 u i u I ( 3 - 45 ) 

2 — 1 

After substituting the estimate R into equa- 
tion (3.44), the final form of the cost function, as 
used in this study, is obtained: 

N 

J{0) = i^(Z i -Y I ) T R _1 (Z J -Y t )+ Constant (3.46) 
1=1 

The cost function given by equation (3.46) is 
the same as that used in an output error technique 
(ref. 12) except the measurement noise covariance 
matrix is used as a weighting matrix. The prob- 
lem is now in the form of an unconstrained optimiza- 
tion problem where the cost function given in equa- 
tion (3.46) must be minimized with respect to the 
unknown parameter vector 0. The unknown param- 
eters determined by this method are, for this study, 
the airplane stability and control derivatives and trim 
coefficients. In addition, measurement noise statis- 
tics (weighting matrix) and parameter standard er- 
rors are determined. 
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The standard errors determined are the Cramer- 
Rao (CR) lower bounds, providing a measure of the 
maximum achievable accuracy for the parameters. 
These are defined by the CR inequality 

E{{0 -0){0- 0) T ) > -E | J (3.47) 


have shown that the ML method produces acceptable 
estimates in many situations. The large sample prop- 
erties of an ML estimate are summarized as follows 
(Cramer, ref. 31, derives these properties): 

1. Asymptotically unbiased: 

lim E{0} = 0 

N^oo 


where 


-E 


d 2 L{0) 

d0d0 T 


= M 


(3.48) 


and M is usually referred to as the Fisher informa- 
tion matrix. It is assumed in this study that the 
approximated Hessian matrix H from the optimiza- 
tion procedure is a good approximation of the Fisher 
information matrix. The solution using a gradient 
optimization scheme generally has the following form 
(ref. 13) for the fcth iteration 


2. Asymptotically efficient: 


E{{0-0){0-0) T } > -E 


d 2 m \ 1 

d0dO T j 


3. Consistent: 


lim Pr{(0 — 0) < e} = 1 

jV— KX> 


with e arbitrarily small 


Ok+1 = h ~ H fc 


(3.49) 


4. Asymptotically normal: 


where 


Sfc 


dJ(0) 

d0 


and Hr;M 


100 


2. Statistical Properties of ML Estimates 

The maximum likelihood method is popular, es- 
pecially in flight test data analysis, because of the 
excellent large sample properties of its estimates. Al- 
though ML estimates do not possess optimal proper- 
ties for small samples, sampling experiments (ref. 30) 


0 = N{0,M~ 1 ) 

Asymptotic unbiasedness and consistency are very 
similar. However, consistency implies that if an esti- 
mator is consistent for 0 , it is also consistent for any 
well-behaved function of 0. Thus, consistency is more 
significant than unbiasedness. Asymptotic efficiency 
is a statement involving the CR inequality; therefore, 
for large samples the CR lower bounds are obtained. 
Asymptotic normality states that as the sample size 
gets very large, the estimates 0 approach a normal 
distribution with mean and covariance given by 0 and 
M —1 , respectively. 
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IV. OPTIMIZATION TECHNIQUES 

The ML parameter estimates are obtained by 
solving an unconstrained, nonlinear optimization 
problem; that is, find 0* which minimizes the cost 
function J{0). The necessary and sufficient condi- 
tions for this problem are as follows: 

1. J{0) is differentiable at 0* . 

2. VJ(0*) = 0. 

3. V 2 J(0*) > 0. 

The theory for solving unconstrained, nonlinear 
optimization problems is often based on the assump- 
tion that the cost function J(0) is a quadratic func- 
tion of 0. This approximation provides a more 
tractable theory and allows basic theorems and prop- 
erties of the optimization methods to be readily es- 
tablished. Corresponding theorems for solving gen- 
eral nonlinear functions of 0 are very difficult to 
prove. However, techniques developed using the 
quadratic assumption are still very effective for non- 
linear functions. Many techniques for solving nonlin- 
ear minimization problems are developed from prac- 
tical experience. 

Optimization techniques for unconstrained prob- 
lems can be divided into two categories: deriva- 
tive methods and search methods. Derivative meth- 
ods may be further classified by the order of the 
derivatives used; search techniques can be divided 
into direct search and random search. Some tech- 
niques combine search and derivative methods; how- 
ever, these hybrid methods are not considered in this 
study. 

The choice between optimization categories de- 
pends on the particular problem. Search methods 
determine the optimization path solely from cost 
function evaluations and, therefore, do not require 
as much algorithm preparation as needed when us- 
ing sensitivity equations. Search methods also do 
not need the regularity and continuity conditions for 
the cost functions that derivative methods need. In 
many unconstrained, nonlinear programming prob- 
lems, however, derivative methods converge faster 
(ref. 13), particularly for estimation problems involv- 
ing dynamic systems. This was demonstrated for air- 
plane systems both in this study and in reference 32. 

Various derivative techniques are available for a 
variety of nonlinear programming problems; how- 
ever, no one technique is best for all problems. For 
example, the steepest descent method works better 
away from the minimum, whereas Newton’s method 
works better near a minimum. A compromise be- 
tween these two techniques is the modified Newton- 
Raphson method (MNR). MNR belongs to a class of 
methods known as quasi-Newton or large step gradi- 
ent methods; these methods approximate the Hessian 


matrix or its inverse while using only first derivative 
information. 

The derivative information can be computed in 
a variety of ways. For dynamic systems, integrat- 
ing analytically derived expressions (sensitivity equa- 
tions) for the derivatives is probably the most ac- 
curate as well as the most time consuming. One 
alternative is a numerical approximation scheme. 
Finite-difference methods are often used because 
they eliminate the additional burden of deriving 
and incorporating sensitivity equations into the al- 
gorithm. However, the finite-difference methods re- 
quire about the same computational effort as inte- 
grating the sensitivity equations. Another option is 
the proposed surface-fitting method of the MNRES 
algorithm presented in section IV.B. 

A. Commonly Used Methods 

Two optimization schemes, representing the two 
main categories of methods, are selected in this study 
primarily to provide a benchmark comparison with 
MNRES. They are commonly used in aircraft es- 
timation and control problems and, therefore, are 
good indicators of the relative merit of MNRES. The 
two optimization methods are the flexible polyhedron 
search (FPS) and the modified Newton-Raphson 
method (MNR). More details are provided in refer- 
ences 13 and 33 for FPS and in references 12 and 
34 for MNR. In a variation of MNR the derivative 
information is computed by using finite differences 
(refs. 16 and 17). Both the finite-difference form and 
the sensitivity equation form of MNR are used in this 
study. 

1. Flexible Polyhedron Search 

Since FPS has been found to be advantageous 
in some aircraft design and control applications 
(ref. 35), it may be a good candidate for reducing 
computational demands in aircraft estimation prob- 
lems. FPS avoids derivative calculations, where the 
quasi-Newton methods spend most of the computa- 
tional time. The algorithm is independent of model 
form and, thus, is readily applicable to any aerody- 
namic model. 

Consider the unconstrained optimization problem 
of minimizing J(0), a scalar function of n p variables. 
The FPS method uses a flexible polyhedron surface 
or simplex with n p + 1 vertices, each defined by a 
vector 0. The vertex 0 #, producing the highest 
value of J(0), is projected through the centroid of the 
remaining vertices to define a new vertex. This new 
vertex, and the remaining ones without 0jj, form a 
new polyhedron. This operation is called a reflection. 
Figure 2 shows two steps in this process for the case 
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ABC Initial simplex 
BCD New simplex 
F Reflection point 
E Contraction point 
G Extension point 


Figure 2. Two steps of a two-dimensional flexible polyhe- 
dron search. 


zero mean and uncorrelated. This leads to a nonlin- 
ear estimation of unknown parameters. Consider the 
system equations (repeating eqs. (2.3) and (2.4)) and 
the measurement equations, 


X s = /(X s ,U,0) X s (0) = X g0 (2.3) 

Y = h(X. s , U, 0) (2.4) 

Zf = Yj + Vj- (i = l,2,...,N) (4.1) 


with two unknown parameters. If the new vertex 
produces a lower cost than 0 £ (the vertex producing 
the smallest J(0)), then an expansion takes place 
and a new vertex is located farther out along the 
same projection. Similarly, if higher costs are found, 
a contraction takes place. The minimum of the 
cost function is found by repeatedly deleting the 
point having the highest value of J(0 ) and adding 
new projected points that produce lower J(0). The 
flexible polyhedron is able to adapt to the shape 
of J{0) by stretching down slopes, contracting near 
minima, and changing direction in curved valleys. 

2. Modified Newton-Raphson Method 

This report is primarily concerned with nonlin- 
ear aircraft estimation problems. Since the MNR 
approach is commonly used for these problems, it 
is included as a benchmark algorithm. Although 
estimating derivatives is computationally burden- 
some, this information enables relatively fast conver- 
gence of the optimization process. In fact, Newton’s 
method converges in one pass for cost functions that 
are quadratic. Hence, Newton and quasi-Newton 
techniques used for estimation problems of dynamic 
systems are expected to converge faster when the 
quadratic approximations for the cost functions are 
valid. Also, these methods provide both step size 
and direction for each iteration. In some problems, 
however, additional control of step size is needed to 
ensure convergence. Since removing the requirement 
of solving sensitivity equations is desirable, the MNR 
algorithm in this report uses a simple finite-difference 
method except when otherwise noted. This is not 
too costly in terms of computational time (refs. 16 
and 17); however, care must be taken to obtain the 
derivatives as accurately as possible. 

The MNR and the MNRES algorithms are the 
derivative methods of interest to this study. As dis- 
cussed in an earlier section, the problem is to mini- 
mize the weighted square of the errors between the 
computed model outputs and the actual measured 
outputs. It is assumed that only the measured out- 
puts are corrupted by noise and that the noise is 


with 


£{vj = 0 and £{v t vj } = R<% (4.2) 

where X s , U, and Y are the state, input, and output 
vectors, respectively; 0 is the unknown parameter 
vector; Z z - and v t - are the measurement vector and 
measurement noise vector, respectively, at f = f,-; R 
is a diagonal measurement noise covariance matrix, 
which under the above assumptions is approximated 
by the covariance matrix of the residuals. Without 
process noise the ML cost function to be minimized 
is given by equation (3.46) from which the added 
constant and multiplicative factor of 1/2 are dropped 
without affecting the solution: 


N 

•/(») = E(Zi-Y f ) T R- 1 (Z i -Yf) 


i=l 


(3.46) 


The matrix R is given by equation (3.45): 


1 N 

R = 


N f 

l—l 


(3.45) 


where 

u i = Z i ~Y i {0 o ) (4.3) 


and 0o is the initial estimate of the unknown pa- 
rameter vector. The MNR method accomplishes the 
minimization by expanding Y, the computed output 
vector, about 0o, the initial unknown parameter vec- 
tor. A Taylor series expansion of Y truncated to first 
order is 


Y(0)=Y(0 o ) + ^ 


AO 


(4.4) 


0o 


where AO = 0— 0q. Then substituting this expression 
into equation (3.46) results in a quadratic approxi- 
mation of J. The increment AO is the unknown. 
Differentiating J with respect to 0 and equating the 
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derivative to zero to find a minimum results in 


% = -H + E GiR-'Gi AO = 0 

2 = 1 

(4.5) 
dyk 


N 


N 


where 


G, = 


[d6 e \ 


(4.6) 


and y k and 0£ are the fcth and f?th elements of the Y 
and 6 vectors, respectively. Solving for Ad gives 


( N 


-1 


N 


AO = 


^GjR-iG, ^GfR-V,- (4.7) 


\*' =1 


1 = 1 


This is often written as 


solution). The trade-offs in choosing a surface and a 
surface-fitting method involve the choice between ac- 
curacy of the sensitivities and computational effort. 


1. Algebraic Solution 

The MNRES algorithm is best described by con- 
sidering the computationally least demanding ap- 
proach of using a linear-surface approximation. Ex- 
panding Y (0) in a first-degree polynomial in 0 for 
each point in time and at n p + 1 different points in 
the np-parameter space gives 


yliiO 3 ) — s k0 + -) 1 - s knp 6 j np (4.9) 


AO = — M -1 


dJ 

dO 


0o 


(4.8) 


emphasizing the Fisher information matrix M and 
gradient terms. For the Arth iteration the estimate 
0/c + 1 is given as 0 k+1 = 0 k + A0 k+l . In this 
study, convergence is achieved when AJ k /J k and 
A0 k /0 k are less than 0.001. The sensitivities Gj are 
determined separately from the above steps. This 
may be done by integrating the sensitivity equations, 
by using a finite-difference approximation, or by 
using MNRES. 


B. MNRES Method 

The MNRES method developed in this report is 
essentially an MNR optimization algorithm with an 
efficient method for estimating sensitivities. As in 
the ML/MNR algorithm previously described, the 
same equations (eqs. (4.1)-(4.8) and (3.45)-(3.46)) 
apply for ML/MNRES; however, the sensitivities Gj 
are computed by using slope information from local 
surface approximations of Y ( 0 ) . The approximations 
are made near the series expansion point of equa- 
tion (4.4). The sensitivities obtained from the fit- 
ted surface are determined with less computational 
effort than that required by either a finite-difference 
method or integration of analytically determined sen- 
sitivity equations. 

The MNRES algorithm is readily optimized for 
a particular application because the user can select 
both the type of surface and the method of fitting the 
surface. Two very practical types of surfaces in aero- 
nautical applications are nth-order polynomials and 
splines. Two efficient methods of fitting the surface 
are by solving n p simultaneous equations for n p un- 
knowns (algebraic solution) and by solving a redun- 
dant set of equations for n p unknowns (least squares 


where i indicates the fth point in time; k indicates 
the fcth element of the output vector Y (0); and j 
indicates one of the n p + 1 sample points used to fit 
equation (4.9) to Y(0). Note that 


yi(0 j ) = y k (0) (4.io) 

at each of the n p + 1 points. The sample points are 
chosen by allowing a small perturbation of each pa- 
rameter around the point where the sensitivities are 
desired. Alternatively, the perturbation size can be 
selected to reflect the relative significance of each pa- 
rameter to the model. This allows for larger pertur- 
bations of the less sensitive parameters and smaller 
perturbations for the very sensitive parameters and 
thus provides higher quality derivative calculations. 
This alternative is discussed further at the end of this 
section. The slopes s k i to s knp are the desired sen- 
sitivities, {dy k /d0(>)i, and s k Q is the value of y k (0) 
at the series expansion point of equation (4.4). Note 
that because this is a linear surface, the slopes are 
constant over the surface and need not be evaluated 
specifically at s k q. If a higher degree polynomial is 
fit to y k {0), the slopes vary across the fitted surface 
and, therefore, must be evaluated specifically at s k 0 . 
Consider the matrix representation of equation (4.9) 
for the first element of Y and for the n p + 1 sample 
points at time tf 


Yh = XS u (4.11) 


Note that y k is the fcth element of the output vec- 
tor Y and y 3 k is the j th element of the surface-fitting 
vector, Y k . Matrix X contains n p + 1 rows defining 
the n p + 1 sample points. Expanding equation (4.11) 
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to show the vector and matrix elements gives 


y (0) 


[yOl 

i'll 


-i 


*8 

y\i 

= 

i 

o\ 

^2 

.yii. 


,i 

6 J* 

fi n 

u 2 



s io’ 


sn 


. s ln . 


(4.12) 


where n = n p for the computationally least demand- 
ing approach, or algebraic solution. In general, n 
may be greater than n p producing a least squares so- 
lution. Since sio is known, equation (4.12) can be 
simplified. The first line in equation (4.12) can be 
eliminated by subtracting it from the other n p equa- 
tions. Thus, 

AY lt - = AXS lt - (4.13) 

or in matrix notation 


'yli-y 0 u 
y\i - y°u 

.yu-y 0 ii- 

where 


A0} 

A 0\ •• 

i-H 

CCS 

<1 

A 0\ •• 


A q •• 

A9J: 

1 

■'-5*5^ 

II 


a*A- 

A*A 


s 11 
5 12 


A^J 


LSlnl i 
(4.14) 


Thus, at time tj, the sensitivities for the first element 
in Y are given by 


S u - = (AX) _1 AY 1 j (4.15) 


Because the AX matrix is independent of time, the 
sensitivities can be calculated rapidly during each 
iteration of the algorithm. This is a key factor in 
reducing the computational effort of the algorithm; in 
effect, the integration of the n s n p (number of states 
times number of parameters) sensitivity equations 
has been replaced by a set of n 0 (number of outputs) 
matrix multiplications. 

Figure 3 shows, geometrically, two iterations for 
the case where 0 is two-dimensional and a linear 
surface is used to fit a scalar y. The expansion at 
time £j is 


y\ - y? 
vi - y,° 


t> 

A 0\' 

'si‘ 

A0j 

A0|. 

. s 2 . 


(4.16) 


During the first iteration, this expansion requires 
that y(0) be evaluated at n p + 1 = 3 points: y°, y 1 , 
and y 2 . Computationally, the first iteration is the 
most costly phase of the MNRES algorithm. Each 



Figure 3. Linear-surface fit for two iterations of MNRES. 

evaluation of Y requires that the equations of mo- 
tion be integrated. The linear surface (indicated by 
the solid triangle in fig. 3) is fit and the slopes (sen- 
sitivities) are thereby determined. The algorithm 
proceeds, as in the ordinary MNR method, by us- 
ing equation (4.8) to obtain 


0 r +i = 0 r + A0 r+ i (4.17) 

The estimated sensitivity values s ^ are used to de- 
fine the elements of matrix G t in equation (4.7). 
The new Y is evaluated (by integration of equations 
of motion) at 0 r + 1 to get y 3 (0). At this point the 
MNRES algorithm has reduced the sensitivity prob- 
lem to solving a set of simultaneous equations. This 
is done by eliminating the 0 3 in X which produced 
the greatest value of J(0 ) and replacing that infor- 
mation with the newest estimate of 0. In the example 
in figure 3, y° was assumed to be the high cost point 
and so was eliminated to obtain the new fitted sur- 
face (indicated by dashed triangle). The slopes of the 
new surface provide the sensitivities for the MNRES 
algorithm to proceed. In this scheme, a check should 
be made to ensure that the new Y J (0 r+ i) produces 
a smaller value of J(0). Sometimes, step-size control 
or complete restarting may be needed. Note that 
initialization of the algorithm requires that n p + 1 
integrations be performed for the n p + 1 trajectories 
Y 3 . After this initializing pass, only one integration 
of the system equations is needed to evaluate the cost 
J ( 0 ) and outputs Y ( 0 ) and to update parameter es- 
timates for each iteration. 

As mentioned previously, in practice it is bene- 
ficial to choose the perturbation size in a different 
fashion from that used in a simple finite-difference 
method. Simply using a 1-percent perturbation on 
each element of 0 to obtain the corresponding per- 
turbation in each element of Y ( 0 ) is not optimum for 
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derivative calculations. Experience has shown that 
it is beneficial to use perturbation sizes that reflect 
the importance of the parameter to the model. By 
computing the sensitivities for each parameter 

and then letting the perturbation sizes be scaled in- 
versely proportional to the normalized ratios of sensi- 
tivities, more accurate derivative information can be 
obtained. Of course, this applies only when an ini- 
tialization, or “restart,” is needed. The fundamental 
issue is that the less sensitive a parameter, the larger 
the perturbation necessary to obtain an appropriate 
size response in the outputs. This approach could 
also be applied to an MNR method. Theoretically, 
the same derivative should be obtained for any suf- 
ficiently small perturbation in 0; however, because 
of the sometimes widely varying sensitivities of the 
parameters as well as the numerical precision limi- 
tations, it is beneficial to vary the perturbation size 
according to the aforementioned rule. The sensitivity 
defined as 6jM^ was introduced in reference 36 and 
used again in reference 37 as a means of quantifying 
the significance of a parameter to the model. 

2. Least Squares Solution 

The least squares approach to fitting the surface 
Y(0) offers another advantage if a recursive least 
squares method is used. The recursive method re- 
duces the storage requirements from n p + 1 sets of 
output time histories to just two time histories: one 
corresponding to the new response predicted by the 
most recent estimate of 0 and the other correspond- 
ing to the outgoing 0 that produced the highest 
cost. The penalty for this advantage is the need 
to integrate equations of motion twice per iteration; 
still, the least squares solution requires substantially 
less computational effort than that required with the 
usual MNR method. 

When using the recursive least squares approach, 
only two changes are made to the MNRES algorithm 
just described. The first change is in the calculation 
of the AX matrix, and the second change is in 
the sensitivity calculation. The development of this 
formulation begins with equation (4.9) in condensed 
form (everything discussed up to this equation in the 
previous development applies here): 

yii = xjs ki ( 4 . 18 ) 

Simplifying the notation by dropping the i subscript 
and writing the matrix form of the equation (which 
removes the j superscript) gives 

Y k = XS k (4.19) 


The least squares solution for the sensitivity vector 
is 

S k = (X r X) ~ 1 X. T Y k (4.20) 

Now, dropping the k subscript and letting the fol- 
lowing apply to the fcth element of the output vector 
Y, a recursive relation can be defined for the r + 1 
iteration 

P r+1 = (X^X,.)- 1 (4.21) 

and the updating equation can be defined as 

P r+l = (P r _1 + a^a r - af a°)~ l (4.22) 

where the row vector a r is the new set of 0 to 
be included and a° is the outgoing set of 0, which 
produced the high cost J. The recursive relation for 
S, which has already been derived in section III.A, is 
now 

S r -)_i = S r — P r+i [a° y° — aj y r 

+ (a^a r - a° T a°)S r ] (4.23) 

where y° and y T are the outgoing and incoming 
elements of Y k , respectively, at the rth iteration. 
With the new sensitivities determined, the algorithm 
proceeds as before. 

3. Properties of MNRES 

Properties of MNRES are discussed in compari- 
son with the commonly used MNR algorithm. Thus, 
convergence characteristics and computational ad- 
vantages and disadvantages of MNRES are compared 
with those of a well-known benchmark. 

Convergence of NR or MNR algorithms, both 
with and without finite-difference derivatives, has 
been well documented (ref. 13). Convergence of 
MNRES can be shown, at least heuristically, by con- 
sidering several details. First, the MNRES method 
is still fundamentally an NR method or, for this 
study, an MNR method. The only critical difference 
is that the derivatives are approximate, which makes 
MNRES closer to MNR with numerically determined 
derivatives. Second, note that fitting a first-degree, 
rip-term polynomial to n p + 1 data points is equiv- 
alent tq a simple finite-difference method. In effect, 
as A0 J (the distance between points on the fitted 
surface for MNRES) becomes small enough, the sen- 
sitivities become identical to those given by a sim- 
ple finite-difference method, regardless of the actual 
functional representation of Y(0). The MNRES al- 
gorithm simply relaxes the accuracy of the sensitiv- 
ities in order to reduce substantially the integration 
requirements; the degree of relaxation varies during 
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the optimization process but can be controlled by 
limiting step size. 

The relaxation of sensitivity accuracy generally 
appears to be a very beneficial trade-off for Newton- 
Raphson algorithms. Before considering this relax- 
ation of sensitivity accuracy, note that during an 
MNRES optimization there are two times that the 
MNRES scheme computes sensitivities, which are 
very close to those computed by a finite-difference 
scheme. These times are, first, during the initial- 
ization, or first pass, of the algorithm and, second, 
toward the end of the optimization process. During 
initialization n p different 0 are chosen (a perturba- 
tion on each element of 0 is sufficient) to give n p 
different response time histories Y (0). The surface 
given by Y ( 0 ) is fitted. The initial 0 ] can be chosen 
such that 

\{0i-O o )/O o |<1 (4.24) 

for each j. For the algebraic solution form of 
MNRES this is completely equivalent to a simple 
finite-difference scheme, and for the least squares 
form it is a very close approximation. In this study, 
the same A was used in the MNR with finite- 
difference derivatives as that used in the initialization 
of MNRES. This was done for comparison purposes; 
in practice, the choice of perturbation size for 9 may 
be very different, as discussed later. Toward the end 
of optimization the MNRES scheme again becomes 
equivalent to a simple finite-difference scheme since 
the AO become very small. This forces the surface 
that is fit to Y(0) to become very small; thus, the 
slope information is computed for a surface fit to a 
very small area. 

The relaxation of sensitivity accuracy occurs be- 
tween the two stages discussed above, that is, after 
initialization and before convergence. During this 
part of the optimization, large AO may occur; this is 
characteristic of NR, MNR, or MNRES algorithms. 
For MNRES, unlike NR or MNR, these large steps 
cause the surface fit area to expand, so that the slopes 
or sensitivities no longer approximate the slope of 
the Y ( 0 ) surface at a “point,” that is, no longer ap- 
proximate the limit requirements of a derivative, but 
rather average the slope over a larger area. This is a 
critical time period for the MNRES optimization. 

Three factors aid in preventing divergence during 
the critical time period. The first factor is that as 
the optimization process advances, MNRES contin- 
ually eliminates values of 0 which are far from 9*, 
the optimal solution. This, in effect, contracts the 
expanding surface that is fitting Y(0) and balances 
the expansion process. As the updated estimates of 
0 become close to 0*, the contraction process domi- 
nates and slope (sensitivity) information approaches 


that given by a finite-difference scheme. The sec- 
ond factor is that Newton’s algorithm and variations 
like NR, MNR, and MNRES advance more quickly 
as the quadratic approximation of the cost function 
improves; moreover, the Newton algorithm converges 
in one step for a quadratic cost function. Since the 
quadratic approximation of the cost function gener- 
ally improves the closer 0 gets to 0*, and since initial 
estimates of 0* are often given by a least squares 
procedure or by a knowledgeable user, Oq tends to 
be “close” to 9*. Thus, for aircraft estimation prob- 
lems, MNRES generally begins in a region conducive 
to convergence. The third factor is that step-size con- 
trol logic can always be incorporated. Carried to the 
extreme, MNRES could always be forced to produce 
the same derivative estimates as a finite-difference 
method. Of course, convergence would be very slow 
because of the very small steps. In practice, one lets 
the algorithm take steps determined by the NR logic 
(as done in this study); then if a convergence problem 
develops, step-size control can be incorporated. 

The computational advantage of MNRES is tied 
to two primary factors. The first factor is the number 
of unknown parameters n p . Both MNR and MNRES 
must integrate n s + n a n p differential equations on 
the first iteration; after that, however, MNRES inte- 
grates only n s state equations during each pass and 
MNR continues to integrate n s state plus n s n p sen- 
sitivity equations. It appears that as n p gets larger, 
so does the advantage for MNRES. A limiting factor 
in MNRES is in equation (4.15), where AX must be 
inverted. This n p x n p matrix becomes more difficult 
to invert as n p gets larger and, unfortunately, is made 
up of very small numbers as the optimization process 
converges. Also, note that the information gained 
during each pass is not equivalent between the two 
methods. MNRES performs less computation dur- 
ing each pass and, consequently, has less information 
with which to update the estimates. However, when 
sufficient passes have been performed to make the 
work done by MNRES equal to MNR, MNRES has 
already performed n p + 1 parameter updates. This 
allows MNRES to step more quickly toward the final 
solution. MNRES achieves convergence faster than 
MNR as the cost function becomes more quadratic. 
The second factor and primary reason for the success 
of MNRES has to do with the degree to which the 
cost function can be approximated by a quadratic 
function. The quadratic approximation is inherent 
to the Newton type of optimization scheme and, 
therefore, both MNR and MNRES improve in per- 
formance as the quality of this approximation im- 
proves. However, convergence occurs more quickly 
with MNRES. This makes sense in light of the way 
convergence takes place. 
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Convergence takes place through an iterative pro- 
cess in which estimates of the unknowns are updated 
during each iteration. The updates are estimated by 
equation (4.8), which is the product of the informa- 
tion matrix and the gradient of the cost function. 
It is well-known that convergence can be speeded up 
by holding the information matrix constant (see, e.g., 
ref. 14) for a limited number of iterations. This elim- 
inates the need to integrate the sensitivity equations 
for a limited number of iterations; note that integrat- 
ing the sensitivity equations accounts for most of the 
computational effort. There are two choices, each 
representing one extreme, for optimization: (1) inte- 
grate the sensitivity equations to obtain the most ac- 
curate derivative information for each iteration (this 
is the most costly in terms of computational effort); 
or (2) hold the information matrix constant for a 
limited number of iterations (this is the least costly 
in terms of computational effort and the least desir- 
able in terms of convergence, since there is no way to 
know for how many iterations the information matrix 


can be held constant without causing divergence). A 
compromise between these extremes is preferred. 

MNR.ES provides this compromise in a very ef- 
ficient manner. A trade-off between computational 
effort and sensitivity accuracy is made automatically 
by MNRES. By using the surface-fitting technique, 
only the state equations need to be integrated dur- 
ing each pass, and this information is incorporated in 
the solution with relatively little computation. The 
sensitivities are only approximated in this process; 
however, their accuracy is controlled sufficiently to 
allow convergence. 

The primary disadvantage of using MNRES comes 
from the computer memory required; n p + 1 sets of 
output variable time histories must be retained. The 
recursive least squares method discussed earlier re- 
duces this storage requirement to two sets of time 
histories; however, the computational effort increases 
from n a to 2 n s equations to be integrated during 
each pass. The user’s computer system would dic- 
tate which approach is more appropriate. 
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V. CONFIDENCE INTERVAL 
ESTIMATION 


(0en) 


(5.2) 


Confidence interval estimation (CIE) is an inte- 
gral part of the parameter estimation problem. Point 
estimates of parameters without any qualifications to 
indicate their accuracy are of little value. An interval 
estimate that incorporates both variance and confi- 
dence level information provides a complete state- 
ment of the estimate quality. Although the Cramer- 
Rao lower bound is commonly used to qualify ML 
parameter estimates, it is well-known that in air- 
craft applications these bounds do not accurately re- 
flect the true parameter variance. They are usually 
too optimistic (ref. 20). The difference between the 
lower bound and the actual parameter variance can 
be due to incorrect assumptions about measurement 
and process noise, bias errors in the estimates, or 
modeling error. However, the nonlinearity of the es- 
timation problem appears to contribute significantly. 
In this chapter a method is discussed for determining 
confidence intervals by analysis of the confidence re- 
gion contours using a search scheme. In addition, a 
measure of nonlinearity is developed to further char- 
acterize the problem. 

A. Confidence Regions 

Confidence regions are described by a surface in 
parameter space representing a certain confidence 
level. The surface is defined by a statistic reflecting 
the distribution of error in 0. From the distribution 
of the statistic, a statement can be made about 
the probability of the statistic being in a certain 
interval I. Assuming that the relationship between 
the statistic and the parameters can be described, a 
further statement can be made that the parameters 
are contained in region Rc with the same probability. 
Region Rc reflects the variation in 6 as the statistic 
varies in interval J. The above procedure is the 
general process by which any confidence interval or 
region is defined. This definition obviously varies 
according to the definition of the statistic. A useful 
statistic for composite hypothesis tests is created 
from the ratio of likelihood functions. 

Let Zi,Z 2 ,...,Zn be N independent random 
variables with probability density functions 


N 

L(n) = Y[f(Zi,0) 

i = 1 
N 

L(oj) = Hf(Z i ,0) (0euj) (5.3) 

i=l 


and denote the maxima of the likelihood functions as 
L(Cl) and L(6j)\ the likelihood ratio is defined as 


A = 


L{U) 


(5.4) 


This ratio forms a statistic having a value between 
0 and 1, since the numerator is limited by the Hq 
hypothesis. The value of A reflects the degree to 
which the Hq hypothesis is accepted, and therefore, 
A can be used as a statistic to test the hypothesis. 
If a probability density function can be defined for A 
and the relationship between A and 0 can be solved, 
then a confidence region R c can be defined. With 
the confidence region determined, the confidence in- 
tervals (extrema of parameters within the confidence 
region) can be defined. 

The confidence region Rc provides an exact de- 
scription of the parameter error bounds. However, 
for the general nonlinear estimation problem, an ap- 
proximation may be involved in defining the confi- 
dence level associated with R c . To resolve this prob- 
lem, Beale (ref. 23) recommended that the statistic 
for the linear estimation problem be used along with 
a correction factor to account for moderate nonlin- 
earity in the model. Since this approach for solv- 
ing the nonlinear problem is based on a correction to 
the linear problem, the following development begins 
with the linear case. 


1. CIE for the Linear Estimation Problem 

The estimation problem is defined to be linear if 
the model equations are linear in the unknown pa- 
rameters; the state, input, and response variables 
may or may not appear linearly in the model equa- 
tions. The form of the linear estimation model (single 
output) is given by equation (3.5), repeated here with 
0 as the vector of unknown parameters (the number 
of measurements is taken to be N for this discussion): 


f(Zi,0) (f = 1,2,..., N) (5.1) 

The testing hypothesis is formulated as Hq: 0 e oj, 
where w is a subset of parameter space fi. Define the 
likelihood functions as 


Y = X0 + e (5.5) 

In this linear regression problem, if Y is N(X0,Ia 2 ) 
and the testing hypotheses considered are 

Hq: 0t = 0 (5-6) 
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and the solution for 0 is 


Hy 0 t ?0 (5.7) 


where the subscript t indicates the true value. It can 
be shown that the likelihood ratio has the form 



A = exp {-^2 [ JW-^)]} 

(5.8) 

where 

J (0) = (Y- X0) t (Y - X0) 

(5.9) 

and 

J(0) = (Y -X0) r (Y -X0) 

(5.10) 

The statistic A can be equivalently replaced by 


H = J(0) - J(0) 

(5.11) 


or in practice by the statistic 

N-np J(0) - J(0) 




J(0) 


(5.12) 


where Fotp is the 1 — Q p point of the F(n p , N — n p ) 
distribution and a p is the confidence level. This is 
possible when the model is correct and J(0) — J(0) 
is independent of J(0 ) (ref. 28). In addition, for the 
linear estimation problem, it is known that (ref. 25): 

1. 0 is an unbiased estimate of 0. 

2. The Cramer-Rao bound is reached. 

The confidence region Rc in parameter space can 
now be given as the set of 0 for which 

J(0) - J(0) < n p s 2 F ap (n p , N - n p ) (5.13) 


where s 2 = J(0)/(N — n p ). Once the data are de- 
termined, J(0) is a function of the n p -dimensional 
parameter space. In parameter space the func- 
tion J(0) — J(0) can be represented by the con- 
tours of a surface. The contours are defined by 
J(0 ) = Constant. Again considering the general lin- 
ear problem (single output) in equation (5.5), the 
cost J (0) can be expanded as 

J(0) = (Y - X0) r (Y - X0) 

= Y r Y - 20 T X T Y + 0 T X T X0 (5.14) 


0 = (X r X) _1 X r Y (5.16) 

The ellipsoidal surface with center at 0 is expressed 
in terms of 0 as 

J(0) -J(0) = 0 T X T X0 - 20 T X T Y 

+ 20 T X T Y - 0 T X T XO (5.17) 

Substituting for XJY from the normal equations 
gives 


J(0) - J(0) = (0 - 0) T X T X(0 - 0) (5.18) 

which defines an ellipsoidal surface in the np- 
dimensional parameter space. For the linear estima- 
tion problem the contours form an ellipsoidal surface 
with a single global minimum. The slopes and orien- 
tation of the contour depend on the model and data; 
in addition, they give an indication of parameter cor- 
relation and conditioning of the problem. If the con- 
tours are greatly elongated (indicating that many 
values of 0 give the same cost), an ill-conditioned 
problem may exist. Inadequate data or possibly over- 
parameterization may be the problem. 

With the relationship in equation (5.13), a con- 
fidence ellipsoid in n p -dimensional parameter space 
with center 0 is defined such that the probability is 
100(1 — a p )% that the true parameter point 0 is con- 
tained within the ellipsoid. This can be expressed by 
substituting equation (5.18) into (5.13): 

(0-0) T X T X(0-0) <n p s 2 F ap (n p ,N-n p ) (5.19) 

The confidence limits are determined by realizing 
that the true value of 0 lies inside the ellipsoid if and 
only if it lies between all points of parallel tangent 
planes to the ellipsoid. Therefore, the true value 
lies between the two tangent planes orthogonal to 
a nonzero vector b if and only if (see ref. 38) 

|b T (0 - 0)| < (b^^b) 1 / 2 (5.20) 

where 

H tt = [n p s 2 F ap (n p , N — n p )\~ l Fl (5.21) 


Differentiating equation (5.14) with respect to 0 and 
setting the derivative to zero, the normal equations 
are obtained 


O = -X r Y + X T X0 (5.15) 


and H is the Hessian matrix of J. Therefore, the 
probability is 1 — a p that for all b, 

|b T 0 - b T 0| < [n p s 2 F ap (n p ,N - ^(b^^b)] 1 / 2 

(5.22) 
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This states that the probability is 1 — a p that for all 

= 1) 2, • . . ) ^p), 

\0 t - 9i\ < \fks\fdii (5.23) 

or, expressed in terms of confidence limits, the prob- 
ability is 1 — a p that simultaneously for all 6{, 


and for the multi-output case, 


J{0) J{0 ) — n p s F Qp (n p , Nn 0 — n p ) 

Nn 0 (n p + 2) 

(■ Nn 0 - n p )n p 


(5.29) 


where s 2 is now given as 


6 — \fkoQ < 6 <0 + \fkag (5-24) 

where 

k = TipFotp^Tip, N Tip) (5.25) 

0Oi = s\fd~i (5.26) 

H -1 = [d i:j } (5.27) 


2. CIE for the Nonlinear Estimation Problem 

The nonlinear estimation problem occurs when 
the unknown parameters appear nonlinearly in the 
model equations. In the nonlinear problem, several 
results change from those found in the linear case 
(ref. 28): 

1. Assuming that e is normally distributed does not 
imply that 0 is normally distributed. 

2. s 2 = J(0)/(N — n p ) is no longer an unbiased 
estimate of <r 2 . 

3. There is no covariance matrix in general. 

4. F-tests and lack of fit tests are not valid in 
general. 

Some results remain true, however: 

1. The sum of squares, J{0), still represents the 
square of the distance from (Zi, Z2, Z3, . . . , Z^) 
to a point in the estimation space. 

2. Minimization of J(0) still corresponds to find- 
ing a point in the estimation space closest to 
{Zi,Z2, . . . , Ztf). 

3. Confidence regions can still be defined; however, 
the confidence level is an approximation. 

In reference 23, Beale recommends using a correc- 
tion factor Nfi as a means to extend the confidence 
level definition of A to moderately nonlinear prob- 
lems. For this case, equation (5.13) becomes 


J{0) - J (0) = n p s 2 F ap (rip, IV — n p ) 


x 



N(n p + 2) ' 

(■ N - n p )n p ^ 


(5.28) 


„2 _ J(fi) 

Nn 0 - n p 


(5.30) 


The correction factor represents a measure of 
nonlinearity (normalized curvature) of the solution 
locus near J{0). The method of computing 
for the multivariable aircraft estimation problem is 
discussed in the next section. 

The confidence contours, defined by equa- 
tion (5.29), cannot be determined analytically as 
done in the linear case since the contours are not nec- 
essarily ellipsoidal. The contours may be very irreg- 
ular and may possibly have several local and global 
minima. Figure 4 shows the construction of a con- 
fidence interval for a one-dimensional problem. The 
solid and dashed curves in the figure represent non- 
linear and linear cases, respectively. In parameter 
space, the dashed curve would form a symmetric el- 
lipsoidal surface, whereas the solid curve would vary 
from this shape, depending on the degree of nonlin- 
earity. The confidence interval for the nonlinear case 
is indicated by 0 m ; n and 0 max ; for the linear case, 
the confidence interval is given by the dashed verti- 
cal lines equidistant from 0. The search algorithm 
used in this study for finding the contour bound- 
aries was presented in reference 26. This method 
tests a series of randomly selected points in parame- 
ter space to determine the position of the confidence 
region. Through many iterations, the limits of the 
region are determined by retaining and updating the 
points on the boundary which maximize or minimize 
the unknown parameters. This search algorithm is 
computationally very demanding, even for problems 
with relatively few parameters. 


B. Nonlinearity Measure for Aircraft 
Applications 

The following is a generalization of Beale’s devel- 
opment (ref. 23) of an intrinsic nonlinearity measure 
and its adaptation to the multivariable problem 
of airplane parameter estimation. This is an empiri- 
cal measure of nonlinearity which, in this study, has 
demonstrated some utility in CIE problems. 

If P(0) represents the estimation space (or solu- 
tion locus) in sample space, then P(0) is the point on 
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Figure 4. Construction of confidence intervals for one- 
dimensional problem. 



and, thus, Q ^ is normalized as 


n P s 2 Q $ 


(5.33) 


The values of this measure may still depend on the 
configuration or orientation of the points P(0 W ) rel- 
ative to P(0), but it should not depend significantly 
on the number of points P{0 W ) chosen or on their 
distances from P(0 ) (if not too large). 

To obtain the intrinsic nonlinearity that Beale 
recommends, must be further restricted: is 

the value of N , ^ when the parameters 0 are chosen 
such that T(0) is always at the foot of a perpendicu- 
lar from P{0) onto a tangent plane at P(0). In other 
words, Nt is the minimum value of N^, if the model 
and the experimental design are fixed (see fig. 5). 

The practical computation of the intrinsic non- 
linearity measure is described as follows. The sen- 
sitivities determined during the estimation process 
are needed in advance of the following calculations. 
According to Beale, an empirical estimate of is 


- n p s 2 <^ 


(5-34) 


Figure 5. Two-dimensional sample space with solution locus 
P(6) and tangent plane T(0 ) at P(0). 


the solution locus closest to the measurement Z and 
0 is the point in parameter space which minimizes 
the cost function. If T{0) is defined as a point on a 
tangent plane at P{0) and W different values are cho- 
sen for the vector 0 near 0 (i.e., 0 W , w = 1,2,..., VF), 
then a crude measure of nonlinearity can be written 
as 

W 

\n0w)-T{0w)\ 2 (5-31) 

w=\ 


A graphical representation of these quantities for a 
two-dimensional sample space is shown in figure 5. 
The nonlinearity measure Q ^ is the sum of squares of 
the distances from the points P(0 W ) to the associated 
points T(0 W ) on the tangent plane at P{0). Clearly, 
QtP depends on the number of points P(0 W ) and 
on their distances from P{0). Beale suggests that 
these distances are proportional to the square of the 
distances of P{0 w ) to P{0). If D is defined as the 
sum of squares of the squared distances, then 

W 

D= £ \P(0 W ) - P(0)\ 4 (5.32) 

ty=l 


where n p is the number of unknown parameters and 
s 2 is the sum of squares of residuals. For the multiple 
output case, s 2 is given by equation (5.30). 

The denominator in equation (5.34), D, can be 
formulated as 


w ( N 

D = £ £[Y;(M - YiWfRT 1 

w=l li=l 


x [Yi(0 w ) - Yj(0)] 


2 


(5.35) 


where i is the index of the N data points (time 
variable). 

By letting T{0) be expressed as a first-order 
Taylor expansion while using the sensitivity infor- 
mation from the estimation, Q ^ can be computed 
as 


W N 

QtI>=H £[y i(M - Y i(0) - gm t 

w=l i= 1 

x R -1 [Y t -(0 w ) - Y { (0) - Gii]> w ] (5.36) 

where 


— 


’ dVk 
idO t \ i 


(4.6) 
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ipw = 0 w -0 (5.37) 

Let 

Wi = Y i(B w ) - Y,-(0) (5.38) 

and rewrite equation (5.36) to obtain 
W N 

Qi>=T, - G i ^) r R- 1 (<5Y i - G itp w ) 

w=l i= 1 

(5.39) 

This is now in the form of a standard least squares 
problem. The problem is to find the value of ^ that 
minimizes Q^\ that is, minimize the distance (see 
fig. 5) given as 

d=\P(0 w )-T(0)\ 2 (5.40) 

Therefore, the value of tp that minimizes Q ^ is $ 
given by W sets of least squares minimizations: 



N 

x J2 G i A ~ l6Sr i (®=1,2 IV) (5.41) 

i=l 


Thus 

W N 

Q<t> = E - G !; $ w ) r R~ 1 (^Yj- - G;$ w ) 

w=l i=i 

(5.42) 

For application to the parameter error bound prob- 
lem for aircraft, the following assumptions are made: 

1. Selecting W = 2 n p is sufficient to adequately 
sample the local surface of J(0) near 0. 

2. Selecting A0 as 

|A0;l = a g. (5.43) 

provides a reasonable distance from 0 to sample 
the surface J(0). The goal is to detect the nonlin- 
earity from the tangent plane without going too 
far into the nonlinear range where the curvature 
(based on sensitivity information at 0) no longer 
applies. 

These assumptions represent a proposal for a 
unified approach to computing N#. With this ap- 
proach, results of various studies can be compared 
and the differences between confidence limits based 
on Cramer-Rao bounds and random search can be 
examined. The error bounds determined by random 
search and can also indicate the effect of experi- 
mental design, especially input form and model error, 
on identifiability. 
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VI. APPLICATION TO SIMULATED 
AND REAL DATA 

The examples considered in this chapter demon- 
strate the ML/MNRES algorithm for estimating 
parameters and the search technique for estimating 
confidence intervals of the parameters. These meth- 
ods are compared with commonly used techniques, 
which provide a benchmark for comparison. The 
commonly used techniques are the ML/MNR algo- 
rithm for parameter estimation and the Cramer-Rao 
(CR) bounds for confidence interval determination. 
ML/MNR is used with both analytically and numer- 
ically determined derivatives. The CR bounds, taken 
from the information matrix, are adjusted to the 95- 
percent confidence level. 

Only dynamical systems or airplane estimation 
problems are used in this study rather than classical 
test problems such as Rosenbrock’s function (ref. 13). 
Using classical optimization problems, which usually 
require very little computational time to evaluate 
the cost function, could lead to different conclusions 
about the algorithms. For aircraft estimation prob- 
lems, the bulk of computer time is spent performing 
integrations of the state and sensitivity equations. 
To prevent any bias in the results due to variations 
in programming efficiency or integration techniques, 
only estimation algorithms using the same integra- 
tion method are compared. 

The performance of the methods used in this 
report is evaluated with the following criteria: 

1. Accuracy of estimates 

2. CPU time to termination 

Termination is obtained when parameter and cost 
function fractional changes are computed to be 
within a specified precision. Both cost function 
change A J/J and parameter change A 6/0 are re- 
quired to be satisfied simultaneously to prevent pre- 
mature termination on a plateau where A J 1 


and A0 is relatively large or on a steep slope where 
A 0 <C 1 and A J is relatively large. 

Besides estimate accuracy and CPU time to ter- 
mination, an additional observation provided is the 
number of “equivalent evaluations.” One unit of this 
measure is the amount of calculation required to in- 
tegrate the system equations and to evaluate time 
histories of the output variables. Each method de- 
scribed in this report requires a different number of 
equivalent evaluations to make one update in the pa- 
rameter estimates. This measure provides an indica- 
tion of how efficiently information gained from sys- 
tem integrations is utilized. System integrations are 
the primary computational burden for any estima- 
tion method applied to dynamical systems. 

The examples in this study use both simulated 
data (examples 1-3) and flight data (examples 4-6). 
Except for the first two examples, the parameters es- 
timated are the nondimensional aircraft stability and 
control derivatives conforming to standard NASA no- 
tation. For examples 1 and 2, a simple linear system 
is integrated with an Euler integration method. Ex- 
amples 3-6 involve the airplane problem and use the 
general equations of motion given by equations (2.5) 
to (2.15). These equations are integrated with a 
fourth-order Runge-Kutta integration scheme. For 
comparison purposes, the same integration step size 
and the same computer (Control Data CYBER 175) 
are used in each example. 

For the airplane examples, the ML/MNRES al- 
gorithm is applied through program MAX. MAX is 
a very modular FORTRAN 5 code with dynamic 
memory. The modular format allows aerodynamic 
models or entire system models to be changed eas- 
ily. The dynamic memory capability adjusts core 
memory automatically to the dimensions required. 
A block diagram of the general computing scheme 
for ML/MNRES is given in figure 6. A flowchart 
for program MAX is given in figure 7 with table I 
defining the elements in figure 7. 
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TABLE I. PRIMARY SUBROUTINES FOR PROGRAM MAX AND DEFINITIONS 
OF FLOWCHART BLOCKS IN FIGURE 7 


Name 

Description 

Subroutines 

AERO 

Computes aerodynamic forces and moments with selected aerodynamic model 

COST 

Computes residuals, fit error, R -1 , and cost 

DIFFEQ 

Computes state derivatives from selected equations of motion 

EST 

Computes updated parameter estimates 

HICOST 

Determines whether new estimates reduce cost and updates storage of outgoing 
and incoming parameters and response time histories 

INT 

Main subroutine for management of model inputs and outputs; computes initial 
conditions and input arrays for RK4 and OUTPUT 

MASTER 

Primary subroutine represented by flowchart; handles initializations, input/output 
operations, and memory management 

OUTPUT 

Computes selected output time histories for cost function and plot routines 

RK4 

Fourth-order Runge-Kutta integration scheme 

SENEST1 

Computes sensitivities using a finite-difference method 

SENEST2 

Computes sensitivities using a selected surface-fitting method 

Decision blocks 

FAIL 

lest whether new estimates reduce cost 

PASSES 

Test for maximum number of allowed passes 

PASS #1 

Test for first pass 

RESTARTS 

Test for maximum number of restarts 

R UPDATES 

Test for maximum number of weighting matrix updates 
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Figure 6. Computing scheme for ML/MNRES. 



Figure 7. Flowchart for program MAX. 











A. Simulated Data Studies 

Simulated data offer three advantages for testing 
a new estimation algorithm. The first and most 
important advantage is that the true values of the 
parameters are known; the second is that the problem 
can be well posed and defined without modeling 
error; and the third is that the degree of complexity 
can be selected. In this study, the first two of 
three simulation examples use a single-input, double- 
output, linear, second-order system: 

X S = AX S + BU X s (0)=0 (6.1) 

Y = X s (6.2) 

These examples are used to demonstrate the relative 
speeds and accuracy of the sample estimation algo- 
rithms and to initially indicate the preferred meth- 
ods. The third simulation uses a nonlinear aircraft 
model and simulates varying levels of measurement 
noise found in real flight data. This example demon- 
strates the accuracy of program MAX. 

1. Example 1 

Example 1 demonstrates and compares the FPS, 
MNR, and MNRES optimization schemes in a sim- 
ple ML estimation problem. The MNR method uses 
a finite-difference method to compute derivatives. 
This satisfies one requirement of this study, which is 
to eliminate the need to derive analytical gradients. 
MNR generally performs with about the same speed 
using either numerically determined derivatives or 
integrated sensitivity equations (ref. 16). MNRES 
uses the same finite-difference method as MNR to 
determine sensitivities during initialization; however, 
MNRES uses the recursive least squares form of the 
algorithm during optimization. The recursive form 
is normally used for conserving memory; however be- 
cause of the small memory requirements to store time 
histories in this example, all time histories are stored 
so that only one integration per pass is needed. The 
purpose of this example is to compare the relative 
performance of each method on a problem involving 
a simple dynamic system. 

Example 1 has six unknown parameters. The 
six unknown parameters in equation (6.1) are the 
four elements of the 2 x 2 system matrix A and two 
elements of the control input matrix B. The input 
form was chosen as 

T1 _ j sin t (0 < t < 2 tt) 

1 0 (t > 2 tt) 

and data points were generated every 0.25 second. 
Process and measurement noise was excluded for this 
example. 


Table II shows the true values, initial values, and 
final estimated values of the six unknown parameters. 
Figure 8 shows the input and response time histories. 
All three algorithms accurately converged to the true 
parameter values using only the first 5 seconds of 
data. The MNR method was 30 times faster than 
FPS and MNRES was twice as fast as MNR, or 
60 times faster than FPS. The number of equivalent 
evaluations had similar ratios, that is 715:28:12 for 
FPS:MNR:MNRES, respectively. 

Table II shows clearly that optimization problems 
having reasonable starting values and involving time- 
consuming cost function evaluations should not be 
solved with direct search schemes, such as FPS. This 
result is supported by an independent study using 
aircraft estimation problems in reference 32. Reason- 
able initial values tend to provide a more quadratic- 
like cost function, for which Newton’s method is most 
effective. If reasonable initial values are not available, 
the FPS may be more attractive. In light of the re- 
sults of example 1, further study concentrated only 
on the gradient methods. 

2. Example 2 

Example 2 is provided to compare the robustness 
of MNRES with that of the commonly used MNR. 
The more common form of MNR with analytically 
derived sensitivity equations is used to prevent any 
deterioration of the algorithm due to approximating 
the sensitivities. The system from example 1 is 
analyzed again except measurement noise is added 
and a pulse input is used to excite the system. Two 
cases are considered, each with different levels of 
measurement noise. The noise is zero mean and 
Gaussian with standard errors of 0.0001 and 0.001 
for cases 1 and 2, respectively. Figure 9 shows time 
histories of the input and response variables for the 
two cases. Table III shows the estimation results. 

In case 1 both methods produce equally degraded 
results; however, MNRES still converged to the same 
precision level more quickly. In case 2, with a 
severe noise level and the information limited to 
3 seconds of data, MNRES was unable to converge. 
The results showed that it was oscillating about 
a solution, unable to find a new parameter vector 
that would produce a lower cost. The MNRES 
used on this example had no special step-size control 
logic. The solution that was obtained, however, was 
as accurate as that obtained by MNR, which did 
converge. 

Meeting convergence requirements does not guar- 
antee accurate results; the error in the estimates 
ranged from 5 percent to 130 percent. MNR had 
both the most accurate and the least accurate esti- 
mate. The importance (sensitivity) of a parameter to 
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the model significantly affects the accuracy of the es- 
timate, particularly under these adverse conditions. 
Based on these examples, it appears that MNRES is 
computationally more efficient than MNR while pro- 
viding the same level of accuracy. 

3. Example 3 

In example 3, the accuracy and robustness of 
ML/MNRES are demonstrated by application to a 
nonlinear aircraft simulation with known measure- 
ment noise levels. In addition, program MAX is val- 
idated. For this example and all other aircraft ex- 
amples, the computationally least demanding form 
of MNRES is used to compute sensitivities. This 
form uses the linear surface fit with equation (4.15) 
instead of the recursive least squares form with equa- 
tion (4.23). The simulation involves a nonlinear lat- 
eral model of a general aviation aircraft. 

Three cases are considered: case 1 without any 
measurement noise, case 2 with a representative noise 
level typical of flight data for the aircraft, and case 3 


with twice the noise level of case 2. The standard 
errors of the simulated measurement noise are shown 
in table IV. In each case the noise is zero mean 
and Gaussian. The simulated data were created by 
a fourth-order Runge-Kutta integration with a step 
size of 0.05 second. Table V shows the terms used in 
the nonlinear aerodynamic model to create the simu- 
lation and the parameter estimates obtained through 
analysis of the simulated data. Time histories are 
provided in figure 10 for the three cases. The con- 
trol inputs were the same for all three cases and are 
shown in figure 11. 

Program MAX was applied using two convergence 
criteria: A J/J < 10 -3 and A0/0 < 10 -3 . As ex- 
pected, the estimates of the less easily identified non- 
linear terms, such as C Uar and C { ap , are more quickly 
corrupted as the noise levels increase; however, the 
estimates are still very reasonable and the time histo- 
ries are accurately predicted. Table V shows that the 
MNRES method can be used effectively in estimating 
parameters for nonlinear aircraft systems. 


TABLE II. PARAMETER ESTIMATES AND COMPUTATION 
TIME FOR FPS, MNR, AND MNRES APPLIED TO A LINEAR 
SYSTEM WITHOUT MEASUREMENT NOISE (EXAMPLE 1) 


Unknown parameters, 

True 

Initial 

Final estimated values 
using method — 

0 

value 

value 

FPS 

MNR 

MNRES 

01 

0 

0.01 


0.89 X 10 -7 


02 

-1.5 

-1.6 

-1.5 

-1.5 

-1.5 

03 

1.0 

1.1 


1.0 


04 

-.5 

-.6 

-.5 

-.5 

-.5 

05 

.2 

.25 

.2 

.2 

.2 

06 

.1 

.15 

.1 

.1 

.1 

Cost 





0.11 x 10 _v 

Equivalent evaluation 



715 

28 

12 

CPU time, a sec . . 



2948 


47 


“Central processing unit time on CYBER 175 computer. 
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TABLE III. PARAMETER ESTIMATES AND COMPUTATION TIME FOR MNR AND MNRES 
APPLIED TO A LINEAR SYSTEM WITH TWO NOISE LEVELS (EXAMPLE 2) 

(a) Case 1 

[Noise level with standard error of 0.0001] 


Unknown parameters, 

e 

True 

Initial 

Final estimated values 
using method — 

value 

value 

MNR 

MNRES 

0i 

0 

0.01 

-0.0675 

-0.0684 

02 

-1.5 

-1.6 

-1.471 

-1.471 

03 

1.0 

1.1 

1.009 

1.010 

04 

-.5 

-.6 

-.449 

-.448 

05 

.2 

.25 

.202 

.202 

06 

.1 

.15 

.098 

.098 

Cost 



0.105 x 10~° 

0.105 x 10~ b 

Equivalent evaluation . . 



42 

12 

CPU time, a sec . . . . 



77.54 

24.08 


(b) Case 2 

[Noise level with standard error of 0.001] 


Unknown parameters, 

True 

Initial 

Final estimated values 
using method — 

0 

value 

value 

MNR 

MNRES 6 

0i 

0 

0.01 

-0.705 

-0.410 

02 

-1.5 

-1.6 

-1.228 

-1.549 

03 

1.0 

1.1 

1.757 

.799 

04 

-.5 

-.6 

.159 

-.238 

05 

.2 

.25 

.210 

.251 

06 

.1 

.15 

.087 

.037 

Cost 



0.104 x 10 -4 

0.122 x 10 -4 

Equivalent evaluation . . 



70 

27 

CPU time, 0 sec . . . . 



134.44 

56.05 


a Central processing unit time on CYBER 175 computer, 
6 MNRES did not converge. 
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TABLE IV. STANDARD ERRORS OF SIMULATED MEASUREMENT 

NOISE (EXAMPLE 3) 


Output variable 

Standard errors 

for — 

Case 1 

Case 2 

Case 3 

/?, rad 

0 

0.010 

0.02 

p, rad/sec 

0 

.010 

.02 

r, rad/sec 

0 

.010 

.02 

<f>, rad 

0 

.005 

.01 

a v , g units 

0 

.005 

.01 


TABLE V. PARAMETER ESTIMATES FROM ML/MNRES APPLIED TO SIMULATED 
NONLINEAR AIRCRAFT SYSTEM (EXAMPLE 3) 


Unknown 

parameters, 

Simulation 

Parameter estimates for — 

e 

values 

Case 1 

Case 2 

Case 3 

Cy,o 

0.13 

0.1299 

0.1298 

0.1295 

Cy 0 

-.411 

-.4136 

-.4261 

-.4401 

C Y P 

-.146 

-.1524 

-.1874 

-.2379 

C Y r 

.63 

.6686 

.6070 

.5412 

C Y Sa 

-.053 

-.0618 

-.0733 

-.0872 

C Y 6r 

.075 

.0794 

.0775 

.0751 

Cl, 0 

0 

.0001 

-.0003 

-.0005 

C b 

-.123 

-.1223 

-.1228 

-.1240 

C b 

-.397 

-.3988 

-.4026 

-.4094 

Cl r 

.257 

.2573 

.2409 

.2239 

C ha 

-.182 

-.1815 

-.1778 

-.1755 

C hr 

.007 

.0067 

.0059 

.00497 

C l ap 

Cfl,0 

2.63 

2.6254 

2.519 

2.4359 

0 

-.00005 

-.00008 

-.0001 

Cnn 

0 

.000003 

.0001 

.0005 

Cn p 

-.15 

-.1488 

-.1524 

-.1558 

Cn r 

-.083 

-.0828 

-.0861 

-.0911 

Cng r 

-.0431 

-.0425 

-.0434 

-.0445 

^notr 

1.7 

1.7343 

1.4419 

1.0118 





1 





Figure 8. Time histories of input and response variables for example 1. 



t, sec 


(a) Case 1. 




t, sec 


(b) Case 2. 


Figure 9. Time histories of input and response variables for example 2. 
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B. Real Flight Data Studies 

In this section three examples are considered. In 
each example the model structure and initial param- 
eter estimates were determined with the modified 
stepwise regression (MSR) program of reference 11. 
For examples 4 and 5, the parameter estimation 
problem is solved by using two ML programs. The 
first is program MAX, which uses the ML/MNRES 
algorithm as described in example 3. The second is 
program MAXLIK, which uses an ML/MNR algo- 
rithm. This MNR algorithm integrates analytically 
derived sensitivity equations to obtain sensitivities. 
MAXLIK is a proven code for aircraft parameter es- 
timation, documented in reference 34. In the last 
example, program MAX is used to compute param- 
eter estimates and CR bounds. These bounds are 
adjusted to the 95-percent confidence level and com- 
pared with those obtained with the search method. 

For comparison purposes, both program MAX 
and program MAXLIK use a fourth-order Runge- 
Kutta integration method with the same integration 
step size (0.05 sec in example 4 and 0.04 sec in ex- 
amples 5 and 6). A convergence criterion is set at 
A J/J = 0.001 for both codes. Program MAX nor- 
mally uses an additional criterion restricting the pa- 
rameter change A 0/0] however, it is removed in these 
examples to ensure that both programs converge for 
the same criterion. Both programs use the same 
bias and scale-factor corrections to the flight data 
for each example. These corrections were determined 
by using a compatibility program developed in refer- 
ence 37. The same initial parameter values are used 
by both MAX and MAXLIK. 

1. Example 4 

Example 4 uses flight data from a general avia- 
tion aircraft operating at an angle of attack of 8°. 
The estimation problem involves the nonlinear lat- 
eral model. Table VI presents a comparison of pa- 
rameter estimates and CR bounds from MAX and 
from MAXLIK. Initial values and sensitivities com- 
puted as 6gM£g are also given. Again, there is rea- 
sonable agreement between the two approaches. CR 
bound estimates tend to be a little higher for pro- 
gram MAX. This is probably due to their sensitivity 
to the derivative information. 

Repeating the calculations with program MAX, 
by allowing the sensitivity ratios to be incorporated 
into the initializing derivative calculations, slightly 
improved the overall speed of the algorithm. This oc- 
curred because only one restart was required during 
the optimization process. More improvement would 


be realized in problems where restarting occurs sev- 
eral times. Time histories of the measured flight data 
and predicted response using the estimated model are 
shown in figure 12. Execution times for example 4 
are 793 seconds for program MAX and 1036 seconds 
for program MAXLIK. 

2. Example 5 

Example 5 uses flight data from an advanced twin 
engine fighter operating at an angle of attack of 6°. 
A nonlinear longitudinal model is used. Table VII 
presents a comparison of parameter estimates and 
their CR bounds from MAXLIK and from MAX 
(case 1). Also shown for each program is the time 
to reach convergence expressed in seconds. Program 
MAX converged very close to the same values as pro- 
gram MAXLIK except processing was done in about 
one-third of the time. This example reflects the effec- 
tiveness of MNRES under fortunate conditions (that 
is, where the cost is well approximated by a quadratic 
function and a moderate number of unknowns (11 pa- 
rameters) are determined). The quadratic nature of 
the cost is indicated by a very small value of N The 
value of Nfj, was 0.003 for this example. The mean 
value estimates of MAX are very good and the CR 
bounds are good but tend to indicate a slightly larger 
error bound than those from MAXLIK. 

Although the Fisher information matrix is up- 
dated during each iteration by both programs, pro- 
gram MAX delays updating the weighting matrix 
R -1 until convergence is achieved. The example is 
solved once more by program MAX and the weight- 
ing matrix is updated twice. These results are shown 
in table VII as case 2 for program MAX. The mean 
values are essentially the same since they are inde- 
pendent of the weighting matrix used, except possi- 
bly through some numerical errors. The standard 
errors are slightly closer to the MAXLIK results, 
and further updating brings only very small improve- 
ments. These estimates were obtained by MAX in 
about 46 percent of the MAXLIK processing time. 

3. Example 6 

The sixth example uses flight data from an ad- 
vanced single engine fighter operating at an angle of 
attack of 4°. This example involves a nonlinear lat- 
eral model. In this example, 95-percent confidence 
intervals are estimated using two approaches. One 
approach is based on the CR bound using program 
MAX and the other on a random search technique. 
The 95-percent confidence intervals determined by 
each approach are presented in table VIII. In this ex- 
ample was found to be 0.02; however, it was set 
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to zero for the interval computations. Even with this 
correction, only a very small increase in interval size 
would be obtained. The confidence intervals deter- 


mined by the search method are significantly larger 
than the corresponding CR estimates and indicate an 
asymmetric confidence interval. 


TABLE VI. PARAMETER ESTIMATES, CR BOUNDS, AND 
SENSITIVITIES FOR PROGRAMS MAX AND MAXLIK APPLIED TO 
REAL DATA FROM NONLINEAR AIRCRAFT SYSTEM (EXAMPLE 4) 


Unknown 

parameters, 

0 

Initial 

value 

Program MAX 

Program MAXLIK 


0 

a 

0 

<7 

0j M u 

C Y,o 

0.036 

0.0061 

0.0006 

0.0213 

0.0005 

0.7533 XlO 3 

Cy p 

-.479 

-.4603 

.0075 

-.4608 

.0067 

.3756 Xl0 5 

C Y P 

-.186 

-.1378 

.0485 

-.0604 

.0439 

.9373 xlO 3 

C Y r 

.522 

.6677 

.0289 

.6209 

.0256 

.1915 xlO 4 

C Y ga 

-.08 

-.0504 

.0166 

-.0375 

.0150 

.8129 XlO 3 

C Y6r 

.083 

.0814 

.0043 

.0763 

.0037 

.9922 xlO 3 

C Y afl 

.45 

.4300 

.0592 

.4512 

.0504 

.9399 XlO 3 

Cl, 0 

0 

.0002 

.00005 

-.0001 

.00005 

.9618 xlO 3 

C b 

-.079 

-.0872 

.0015 

-.08 

.0013 

A493 xlO 6 

C b 

-.47 

-.5320 

.0102 

-.4823 

.0085 

.4529 xlO 7 

Cl r 

.187 

.1700 

.0043 

.1543 

.0045 

.6838 xlO 4 

C ha 

-.19 

-.2035 

.0036 

-.1852 

.0031 

.7969 xlO 5 

C hr 

.01 

.00055 

.00024 

-.0012 

.00072 

.1340 xlO 4 

^ l a0 

-.26 

-.2707 

.0116 

-.2105 

.0091 

.6327 xlO 4 

C'n,o 

0 

-.00063 

.00003 

-.002 

.00002 

.1909 XlO 4 

Cn 0 

.04 

.0323 

.00045 

.0329 

.0004 

.1515 xlO 6 

C?i p 

-.056 

-.1043 

.0026 

-.0916 

.0022 

.1400 xlO 6 

Cn r 

-.15 

-.1462 

.002 

-.1534 

.0017 

.5039 xlO 5 

Cns a 

0 

-.0044 

.001 

-.0037 

.0009 

.1780 xlO 4 

Cn 6r 

-.053 

-.0550 

.0003 

-.0532 

.0003 

.9048 xlO 7 


-.39 

-.39 


-.39 


° n (5 3 

.08 

.08 


.08 




“Parameter held fixed. 
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TABLE VII. PARAMETER ESTIMATES, THEIR STANDARD 
ERRORS, AND TIME TO REACH CONVERGENCE FOR 
PROGRAMS MAX AND MAXLIK APPLIED TO REAL DATA 
FROM NONLINEAR AIRCRAFT SYSTEM (EXAMPLE 5) 


Unknown 

parameters, 

0 

MAXLIK 

MAX 

Case 1“ 

Case 2 b 

0 

c er 

0 

c~ 

c tr 

0 

C cr 

C X,o 

0.17017 

0.208 X 10 -3 

0.17078 

0.249 X 10~ 3 

0.17070 

0.212 X 10~ 3 

C X a 

.9464 

.616 x 10 -2 

.9245 

.703 x 10 -2 

.92511 

.555 X lO -2 

°x Se 

.2789 

.404 x 10 -2 

.2754 

.116 X 10 _1 

.2755 

.876 X 10 -2 

Cz,o 

-.83946 

.896 x 10~ 3 

-.84335 

.106 x 10~ 2 

-.84252 

.103 X 10~ 2 

C Z a 

-5.311 

.230 X 10 _1 

-5.197 

.223 X nr 1 

-5.194 

.215 X 10 _1 

CZq 

-18.7 

.162 x 10 

-20.3 

.189 x 10 

-20.5 

.177 X 10 

°z 6e 

-.618 

.264 x 10“ 1 

-.566 

.350 x lO -1 

-.570 

.324 X nr 1 

Cm,o 

-.001251 

.828 X 10“ 4 

-.001610 

.942 x 10 -4 

-.001555 

.937 X 10~ 4 

Cm a 

-.5129 

.102 X 10~ 2 

-.5186 

.115 x 10 -2 

-.5183 

.110 X 10 -2 

Cmq 

-16.95 

.157 

-19.06 

.144 

-18.83 

.144 


-1.3409 

.576 X 10 -2 

-1.4150 

.447 x 10 -2 

-1.4075 

.461 X 10“ 2 

CPU time, d sec . . . 

342 

105 

157 


“Case 1 is for one update of weighting matrix R -1 . 
b Case 2 is for two updates of weighting matrix R _1 . 
c Cramer-Rao bound. 

^Central processing unit time on CYBER 175 computer. 


TABLE VIII. PARAMETER ESTIMATES AND CR BOUNDS FROM 
PROGRAM MAX AND CONFIDENCE LIMITS FROM A 
RANDOM SEARCH TECHNIQUE FOR REAL DATA 
FROM NONLINEAR AIRCRAFT SYSTEM (EXAMPLE 6) 


Unknown 

parameters, 

0 

0 

95% confidence interval 

Cramer-Rao 

bound 

Random search 

Upper 

bound 

Lower 

bound 

C Yp 

-0.77 

±0.27 

0.92 

-0.89 

C Y Sr 

.18 

±.27 

.87 

-.89 


-.228 

±.021 

.042 

-.10 

C b 

-.88 

±.13 

.30 

-.99 

Clr 

-4.20 

±.98 

2.1 

-5.5 

C ka 

-.152 

±.020 

.036 

-.11 

Cn 0 

.0826 

±.0040 

.011 

-.013 

Crip 

-.078 

±.019 

.050 

-.10 

Cn r 

-1.24 

±.19 

.49 

-.76 

Cng r 

-.0860 

±.0064 

.020 

-.021 
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Figure 12. Measured and predicted responses for lateral maneuvering using real flight data (example 4). 
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Figure 12. Concluded. 
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C. Discussion of Results 

The general experience with ML/MNRES and 
the examples chosen for this study indicate that 
ML/MNRES performs better than ML/MNR for 
estimation problems involving dynamic systems such 
as aircraft systems. In general, MNRES should 
perform well in any problems for which the Newton- 
Raphson family of methods is appropriate (that is, 
where the cost can be reasonably approximated by 
a quadratic function). The results of this study 
also indicate that a search technique is needed to 
accurately assess the parameter error bounds in the 
nonlinear estimation problem; for the linear problem 
or problems with very little nonlinearity, the CR 
bounds should agree with values determined by the 
search technique. 

1. Parameter Estimation 

Two goals in this study were to develop an es- 
timation algorithm which, first, eliminated the re- 
quirement to reformulate the algorithm for each new 
model and, second, provided a more efficient method 
of dealing with computationally more burdensome 
nonlinear problems. The first goal has been sur- 
passed through the ML/MNRES algorithm coded in 
program MAX in two respects: (1) it does not re- 
quire the derivation of sensitivity equations to com- 
plete the formulation of the algorithm; and (2) the 
modular form of MAX allows easy application to any 
system. The second goal is achieved in three ways: 
(1) the algorithm provides the user many computa- 
tionally efficient options for approximating the sensi- 
tivities (allowing for variation in accuracy and order 
of derivatives); (2) the algorithm allows for substan- 
tial reduction in memory requirements with only a 
small cost in additional computation; and (3) all the 
above options are readily incorporated because of the 
modular format of MAX. The second goal has been 
demonstrated in the examples and the following dis- 
cussion clarifies the conditions under which this goal 
has been met. 

The first two examples use a simulated linear 
system with and without noise. This system is read- 
ily identifiable except when severe noise and nonop- 
timal inputs are included. Because the initial val- 
ues were relatively close to the final solution, a good 
quadratic approximation of the cost function was 
possible and thus provided a condition conducive to 
convergence. The two Newton-Raphson methods, 
MNR and MNRES, were substantially faster than 
the search method as expected under these condi- 
tions. MNRES, however, capitalized more efficiently 
than MNR on the information obtained from each in- 
tegration of the system equations. Each integration 


of these equations provides information that is imme- 
diately incorporated into the numerical process when 
using MNRES. When using MNR, n p + 1 system inte- 
grations (equivalent evaluations) are required before 
each updating operation; for example 1, the ratio of 
equivalent evaluations was 28:12. The results indi- 
cated both MNR and MNRES to be very fast relative 
to the search technique; thus search methods were 
eliminated from any further study. The MNRES ap- 
proach was a little more than twice as fast as MNR. 

The third example demonstrated that the 
ML/MNRES algorithm was a viable method for a 
nonlinear aircraft estimation problem with both re- 
alistic and heavy noise levels. This example provides 
confidence in the ML/MNRES approach. However, 
since it is a simulated example, it cannot be accepted 
as conclusive. Simulations always provide optimal 
conditons for estimation algorithms because prob- 
lems such as modeling error, bias errors, unknown 
noise spectra, and general data incompatibility are 
not present. 

Unlike simulated data, real flight data often 
present problems (as just mentioned) for any esti- 
mation method; these problems may slow the con- 
vergence process or even stop it. The last three 
examples used real flight data and were specially se- 
lected to reflect differing levels of difficulty for the 
estimation algorithms. Examples 4 and 6 relative to 
example 5 demonstrate a representative range in the 
degree of difficulty (measured by computational ef- 
fort) for the algorithms, and it is no surprise that the 
degree of nonlinearity also varies widely (N^ differed 
by an order of magnitude between examples 5 and 
6, 6 being more nonlinear). ML/MNRES was again 
faster than the benchmark algorithm ML/MNR. For 
the more nonlinear examples, convergence time for 
ML/MNRES was 70 to 80 percent of the time re- 
quired for ML/MNR; in the less difficult problem, 
ML/MNRES required only 46 percent of the bench- 
mark time. These examples give some credence to 
the superiority of ML/MNRES. 

However, these examples also indicate a large 
variability in the superiority of ML/MNRES. As the 
degree of nonlinearity increases, the two methods ap- 
pear to approach the same speed of convergence. The 
computational advantage of ML/MNRES tends to 
be reduced as the nonlinearity increases. A moder- 
ate number of unknowns are used in both example 5 
and example 6 so the advantage due to the difference 
in n s + n s n p integrations per pass and n s integra- 
tions per pass is probably a small factor (see sec- 
tion IV.B). The main factor, however, is the quality 
of the quadratic approximation of the cost function 
which, of course, is directly related to the degree of 
nonlinearity of the cost. Both methods are slowed as 
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Figure 13. Convergence criterion versus CPU time for 

programs MAX and MAXLIK. 

the nonlinearity increases, or the quadratic approx- 
imation becomes poorer, but the MNRES method 
depends much more on the quality of the quadratic 
approximation, since in effect it is an approximation 
of the MNR method. As the nonlinearity increases, 
MNRES loses its advantage over MNR. 

Figure 13 offers a graphical view of the perfor- 
mance of MNR and MNRES in example 5, where the 
quadratic approximation is fairly good. The graph 
shows the value of the convergence criterion versus 
CPU time. Program MAXLIK, using MNR, follows 
a typical convergence pattern; the small oscillations 
before convergence are due to the updating of the 
weighting matrix R during each pass after the cri- 
terion falls below 0.01 in value. This approach with 
MNR has been found to be beneficial in these prob- 
lems. Program MAX, on the other hand, updates 
the information matrix during each pass but holds 
the weighting matrix constant until the first conver- 
gence is achieved. At this point the final param- 
eter estimates are obtained; however, the Cramer- 
Rao bounds (determined from the information ma- 
trix) are not converged. The information matrix M 
does not converge until the weighting matrix R -1 
is updated sufficiently. This is understandable since 
the parameter estimates are asymptotically indepen- 
dent of the weighting matrix, whereas the CR bounds 
depend on the information matrix, which in turn de- 
pends on the value of R. Therefore, two more cycles 
are made to convergence to ensure that the weighting 
matrix has stabilized. Note that the first step in each 
method takes the same amount of time and achieves 
the same reduction in cost; this is expected since ini- 
tializing MNRES requires the same computations as 
the first pass in MNR. 

A key feature of ML/MNRES is the method of 
updating the information matrix. Although both 
ML/MNRES and ML/MNR update the information 
matrix ’uring each pass, they do it quite differently. 


One well-known way to improve the speed of tech- 
niques such as ML/MNR involving the Hessian ma- 
trix (or approximations to it) is to hold the informa- 
tion matrix M constant for one or more iterations 
(ref. 14). This reduces the amount of integration re- 
quired per pass to the same as required in MNRES 
because no sensitivity equations are integrated. The 
number of iterations that M can be held constant is 
unknown and unknowable in advance. So there are 
two alternatives generally known and used. One is to 
integrate the sensitivity equations during each pass, 
requiring maximum computational effort but giving 
maximum accuracy to the optimization process. The 
other alternative is to integrate the sensitivity equa- 
tions during a necessarily infrequent number of iter- 
ations to hopefully increase convergence speed. In 
methods like the ones considered in this study where 
large steps in the optimization process may occur, 
not updating the information matrix is dangerous. 
ML/MNRES offers an effective compromise between 
these two unsatisfactory alternatives. 

The compromise is achieved by updating the in- 
formation matrix during each pass, but incorporat- 
ing only the information obtained from integrating 
the equations of motion once. Thus, updating is oc- 
curring at minimal computational cost. Because of 
the limited information to update the information 
matrix, a suboptimal, but computationally more ef- 
ficient, path is followed to convergence. The result 
is that ML/MNRES requires many more passes to 
reach convergence, but each pass requires much less 
computational effort than each pass in ML/MNR. 
The net result is much faster convergence, depend- 
ing on the degree of nonlinearity of the cost and the 
quality of the quadratic approximation used by the 
method. 

2. Confidence Interval Estimation 

Confidence intervals obtained in example 5 were 
found to be very close to the CR bounds adjusted 
to the 95-percent confidence level. In addition, the 
value of Nfj, was very small and convergence occurred 
relatively quickly. This indicates that the cost func- 
tion was very well approximated by a quadratic func- 
tion. In analogy to figure 4, the construction of 
the confidence intervals for example 5 would place 
the dashed and solid lines virtually on top of each 
other at the error level selected. The quality of the 
quadratic approximation is confirmed by the very 
fast convergence of MNRES relative to MNR. 

Confidence intervals obtained in example 6 were 
found to be 3 to 8 times the size of the CR bound ad- 
justed to the 95-percent confidence level (table VIII). 
This is in agreement with references 20 and 25 on 
the values generally found in analyzing actual flight 
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data. Reference 25 used the search technique and ref- 
erence 20 estimated the confidence intervals by com- 
puting an ensemble average. This also indicates that 
this cost function is much more nonlinear than that 
for the first example and is not well approximated 
by a quadratic function at the error level considered. 
This is confirmed by the relative speed of convergence 
of MNRES and MNR. MNRES, with two updates of 
R, required 85 percent of the time MNR required. 

It appears from the examples considered that a 
primary factor in determining the confidence interval 
size for airplane stability and control derivatives is 
the degree of nonlinearity of the cost function. Other 
factors, such as bias errors, modeling error, noise 
level, and noise spectra, may contribute directly to 
confidence interval size or may manifest themselves 
as part of the nonlinearity of the cost function. 
In this study, the other factors were not tested to 
determine their contribution. 


At present, the random search technique is the 
only method to determine confidence bounds accu- 
rately. Clearly, the CR bounds, which are tied to 
the quadratic approximation inherent in the infor- 
mation matrix, will always differ from the parameter 
variance. This difference will mainly depend on the 
quality of the quadratic approximation. The only 
disadvantage of the search technique is its relatively 
poor convergence rate combined with the large com- 
putational effort required in this type of problem. 
Although Beale’s measure of nonlinearity, N^, was 
designed to correct the confidence level in the CIE 
problem, there seems to be more utility in consider- 
ing N , ^ (or some similar measure) as a way to dis- 
cern whether the lengthy computations of the ran- 
dom search are needed. If very little nonlinearity 
exists, the user can be reasonably confident that the 
CR bounds provide accurate error bound information. 


43 



VII. CONCLUDING REMARKS 

The primary contribution in this study is a 
methodology for solving the nonlinear airplane iden- 
tification problem. The use of a modified stepwise 
regression in conjunction with several testing criteria 
is suggested to determine the airplane aerodynamic 
model structure. This very efficient scheme readily 
accommodates the widely varying model structure 
in nonlinear flight regimes. A maximum likelihood 
scheme with the optimization algorithm developed 
in this study (ML/MNRES) is then recommended to 
obtain optimal parameter estimates. This method 
is more efficient than other commonly used tech- 
niques in airplane estimation problems and provides 
some practical computing options. Finally, a random 
search procedure is required to determine parameter 
confidence limits for nonlinear problems. This is used 
in conjunction with Beale’s measure of nonlinearity 
(adapted to the airplane problem) to make an empir- 
ical correction to the confidence level. It is also used 
to determine whether the extensive calculations of 
the random search are needed to estimate confidence 
limits. 

The new optimization algorithm, MNRES, has 
three advantages over other comonly used tech- 
niques. The first advantage is that the algorithm 
removes the need to derive sensitivity equations for 
each new model; this eliminates the computational 
burden of integrating the sensitivity equations dur- 
ing each iteration of the algorithm and also provides 
much flexibility, allowing the model equations to be 
in any convenient format, such as splines, polynomi- 
als, or a nonanalytic form. Also the quickly varying 
model structure sometimes found in the nonlinear 
regimes is readily handled. The second advantage is 
that the algorithm is effective for a variety of methods 
for fitting the output vector to a surface in parame- 
ter space (needed for sensitivity estimation), the user 
can choose a surface-fitting method best suited to the 
problem. Also storage requirements can be reduced 
with little additional computation. The third advan- 
tage of the algorithm is that it requires less com- 
putational effort than the commonly used modified 
Newton-Raphson (MNR) method. For small prob- 
lems (fewer than 15 parameters to be estimated), the 
reduction can be substantial. For larger nonlinear 
problems, the reduction may be more modest; how- 
ever, improvements may still be significant if data 
quality, signal compatibility, and sensitivity calcula- 
tions are accurate. Even though the application of in- 
terest for this study was an aircraft operating in non- 


linear flight regimes, the approach should be effective 
for many other nonlinear, dynamic systems. Based 
on this study, the ML/MNRES algorithm generally 
performs better and offers more versatility than the 
commonly used ML/MNR algorithm. 

The suggested methodology includes a random 
search technique to obtain parameter confidence lim- 
its for the maximum likelihood estimates. Since the 
nonlinear problem does not lend itself to an explicit 
analytical solution, the search uses a random sam- 
pling algorithm to find the confidence limits; unfor- 
tunately, this method is computationally demanding, 
particularly for problems with a large number of un- 
knowns. Unless sufficient repeated measurements are 
available, it is the only method to accurately deter- 
mine confidence region boundaries in the nonlinear 
problem. Beale’s measure of nonlinearity is used 
to provide an empirical correction to the confidence 
level used by the search. Although this was Beale’s 
intended use, it has little affect on the confidence lim- 
its for airplane applications. However, it was shown 
that the degree of nonlinearity is closely related to 
the degree to which the Cramer-Rao bounds and the 
random search confidence limits agree. Therefore, 
it is recommended that this or some similar mea- 
sure be used to determine the necessity of the search 
calculations. 

If further studies are made with MNRES, it 
should prove beneficial to use more efficient inversion 
schemes than the standard Gaussian elimination used 
in this study. This may improve the algorithm for 
larger numbers of unknowns. Also, further consider- 
ation should be given to defining how the confidence 
intervals and the nonlinearity of the cost function re- 
late to other factors such as bias errors, modeling 
error, input form, and noise spectra. In addition, 
measures of nonlinearity and the best schemes for 
computing them need more investigation. Nonlin- 
earity measures may be useful for reflecting the qual- 
ity of the experiment, since parameter error bounds 
will vary with model error and optimality of the in- 
put form. Finally, significant computational savings 
would be achieved if the confidence limits for the non- 
linear estimation problem could be determined with 
gradient techniques rather than the computationally 
demanding search scheme used in this study. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
October 17, 1985 
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